PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Biomed Eng. Author manuscript; available in PMC Nov 11, 2013.
Published in final edited form as:
PMCID: PMC3822781
NIHMSID: NIHMS506198
Using Point Process Models to Compare Neural Spiking Activity in the Subthalamic Nucleus of Parkinson’s Patients and a Healthy Primate
Sridevi V. Sarma, Member, IEEE,corresponding author Uri T. Eden, Ming L. Cheng, Ziv M. Williams, Rollin Hu, Emad Eskandar, and Emery N. Brown, Fellow, IEEE
Sridevi V. Sarma, Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA;
corresponding authorCorresponding author.
Sridevi V. Sarma: sree/at/mit.edu
She is now with the Department of Biomedical Engineering, Institute for Computational Medicine, Johns Hopkins University, Baltimore, MD 21218 USA
Placement of deep brain stimulating electrodes in the subthalamic nucleus (STN) to treat Parkinson’s disease (PD) also allows the recording of single neuron spiking activity. Analyses of these unique data offer an important opportunity to better understand the pathophysiology of PD. Despite the point process nature of PD neural spiking activity, point processmethods are rarely used to analyze these recordings. We develop a point process representation of PD neural spiking activity using a generalized linear model to describe long- and short-term temporal dependencies in the spiking activity of 28 STN neurons from seven PD patients and 35 neurons from one healthy primate (surrogate control) recorded, while the subjects executed a directed-hand movement task. We used the point process model to characterize each neuron’s bursting, oscillatory, and directional tuning properties during key periods in the task trial. Relative to the control neurons, the PD neurons showed increased bursting, increased 10–30 Hz oscillations, and increased fluctuations in directional tuning. These features, which traditional methods failed to capture accurately, were efficiently summarized in a single model in the point process analysis of each neuron. The point process framework suggests a useful approach for developing quantitative neural correlates that may be related directly to the movement and behavioral disorders characteristic of PD.
Index Terms: Deep brain stimulation (DBS), Parkinson’s disease (PD), point processes, spike trains
The use of chronic deep brain stimulation (DBS) is now a well-established therapy for Parkinson’s disease (PD) [5], [23]. Prior to the placement of the stimulating electrodes in the commonly targeted subthalamic nucleus (STN), microelectrode recordings are routinely performed to help to locate the STN and to insure proper placement of the electrode. The need to record neural spiking activity to insure proper electrode placement gives neurophysiologists a unique opportunity to learn directly about the pathological properties of STN neurons in human patients with PD [1], [4], [5], [10], [26], [29].
To date, the analyses of these unique data typically use several different statistical techniques to characterize the spiking properties of the STN neurons [14]. Short-term history dependence within a spike train is analyzed with interspike interval histograms [37]. Long-term history dependence related to neural oscillations is often studied in the frequency domain using power spectra by transforming the spike train into a continuous-valued signal before computing its Fourier transform [22], [28]. Movement-related properties are determined from a tuning curve computed as the average spike rates across multiple trials in each movement direction [20], [21], [33].
More recently, point process methods have been used to analyze the spike train activity for a broad range of neural systems [6], [25], [30], [35], [42], [43]. Despite the point process nature of STN spiking activity, point process methods are not routinely used to analyze these recordings. Some studies have applied Poisson models to neuronal spiking data of PD patients, but these models fail to capture any temporal dependencies that exist in the spiking activity, such as refractoriness, bursting, and oscillations [26], [29]. Zelnikera et al. [40] analyzed the stochastic structure of the short- and long-term dependencies as two Poisson processes with different rate constants. Although informative, this approach restricted the summary of the temporal dependencies in the data to the differences between the two rate constants. Eden et al. recently demonstrated that a single point process model could be used to characterize the spiking activity in STN neurons recorded from PD patients executing a directed-hand movement task [18]. The model used a generalized linear model (GLM) [27] to represent the point process conditional intensity function (CIF) in terms of both short- and long-term history dependence [35]. The single model captured oscillations, bursting, directional tuning, and thus, obviated the need for multiple different analyses. In addition, the model identified a previously undescribed period of decreased spiking propensity 20–30 ms following a spike immediately prior to movement [18].
To use this recent work to understand more clearly the extent to which the features of the neural spiking activity observed in the PD patient’s are signatures of this pathological condition, it would be ideal to compare these features with similar features derived from spiking activity recorded in healthy human subjects executing the same movement task. Because this is not possible for obvious ethical reasons, we decided to compare the STN neural spiking activity of seven PD patients with STN activity recorded from a healthy nonhuman primate, rhesus monkey, executing the movement task. The nonhuman primate serves as a surrogate control. We use GLMs to formulate point process models [15], [34] to characterize the relative contribution of intrinsic factors (e.g., short- and long-term history effects) and extrinsic factors (e.g., the impact of movement direction) on the probability that the neuron will spike at any given time [35]. Once adequate goodness-of-fit is established, we use the model parameters to analyze the relative importance of bursting, oscillations, directional tuning in characterizing differences in spiking propensity of neurons in the two groups.
A. Human Subjects
Seven patients undergoing deep brain stimulator placement for the treatment of PD were included in the study. All patients had idiopathic PD with a Hoehn–Yahr score [41] of three or higher and had a documented response to L-dopa replacement therapy. All patients received a thorough preoperative neurological exam. Exclusion criteria for surgery included those patients with Parkinson “plus” syndromes, cognitive impairment, active psychiatric disorders, or anatomic abnormalities on MRI [4], [5]. None of the patients had undergone prior surgery for the treatment of PD. Informed consent for the study was obtained in strict accordance with a protocol approved by the Institutional Review Board and the multidisciplinary movement disorders assessment committee at the Massachusetts General Hospital. The decision to offer surgery was based on clinical indications alone, and bear no relation to the patients’ participation in this study. To ensure that the patients were comfortable with performing the behavioral joystick task, they practiced it prior to surgery. Subjects were able to remove their hand from the joystick or stop the task at any time. At all time points before and during surgery, the patients had the clear understanding that their participation was not related to the surgical outcome, and that they could withdraw from the study at any time.
B. Electrophysiology
Anti-Parkinsonian medications were withheld the night before surgery. No sedatives were given prior to or during performance of recordings. A local anesthetic was used prior to the incision and burr hole placement. The stereotactic localization using preoperative MRI and computerized tomography, as well as general techniques of intraoperative microelectrode recordings have been described previously [1], [4], [5]. Single-unit recordings were made from the dorsal-lateral motor subterritory of the STN based on stereotactic localization and reconstructions of the electrode trajectories [1]. The STN has characteristic high- firing rates in comparison to the surrounding structures [2] and has clear dorsal and ventral borders that are evident when reconstructing neuronal activity along the electrode trajectories. Once within the STN, no attempt was made to explicitly select cells based on presence or absence of movement-related activity, or on whether the cells responded to passive and/or volitional movement. This was done specifically to limit the potential for a sampling bias. We used an array of three tungsten microelectrodes, separated by 2 mm and placed in a parasagittal orientation. The electrodes were advanced simultaneously in 50 µm increments using a motorized microdrive (Alpha Omega, Nazareth, Israel). The behavioral paradigm was controlled by a Macintosh G4 computer using custom-made software. Neuronal activity was bandpass filtered (300 Hz–6 kHz) and sampled at 20 kHz. Spikes were sorted offline using a standardized template-matching algorithm (Cambridge Electronics Design, Cambridge, England).
C. Primate Subject
One adult male rhesus monkey (macaca mulatta), a.k.a. “Bohr,” was included in this study. A titanium head post, plastic recording chamber, and scleral search coil were surgically implanted in accordance with guidelines set by the animal review committee at Massachusetts General Hospital. Neuronal activity was amplified, bandpass filtered between 200 Hz–5 kHz, and sampled at 20 kHz. Spikes were stored and sorted offline using a template-matching algorithm (Spike 2, Cambridge Electronics Design). Eye position and joystick position were each sampled and recorded at 1 kHz. Bohr performed 868 trials of a directed movement task described below and recordings were taken from 96 neurons.
D. Behavioral Task
Once the microelectrodes were in the STN, the subjects viewed a computer monitor and performed a behavioral task by moving a joystick with the contralateral hand. The joystick was mounted such that movements were in a horizontal orientation with the elbow flexed at approximately 45°. The behavioral task began with the presentation of a small central fixation point. After a 500 ms delay, four small gray targets appeared arrayed in a circular fashion around the fixation point (up, right, down, and left). After a 500–1000 ms delay, a randomly selected target turned green [target cue (TC)] to indicate where the subject is to move. Then, after another 500–1000 ms delay, the central fixation point turned green [go cue (GC)], cueing the subject to move. At this point, the subject used the joystick to guide a cursor from the center of the monitor toward the green target. Once the target was reached, either a juice reward was given (in the primate case) or a tone sounded indicating that the subject had successfully completed the task (human case), and the stimuli were erased. Subjects were required to return the joystick to the center position before the next trial started. A schematic representation of a single trial is shown in Fig. 1.
Fig. 1
Fig. 1
Schematic of behavioral task trial.
Table I breaks down the number of trials and neurons used in our analysis for each PD subject and for our primate subject.
TABLE I
TABLE I
Distribution of Trials and Recorded Neurons Per Patient and Primate
E. Point Process Model of STN Dynamics
We formulate a point process model to relate the spiking propensity of each STN neuron to factors associated with movement direction and features of the neuron’s spiking history. We use the model parameters to analyze oscillations, bursting, and directional tuning modulations across the entire trial and make comparisons between two subject groups.
A point process is a series of 0–1 random events that occur in continuous time. For a neural spike train, the 1-s are individual spike times and the 0-s are the times at which no spikes occur. To define a point process model of neural spiking activity, in this analysis, we consider an observation interval (0, T], and let N(t) be the number of spikes counted in interval (0, t] for t [set membership] (0, T]. A point process model of a neural spike train can be completely characterized by its CIF, λ(t|Ht) defined as follows:
equation M1
(1)
where Ht denotes the history of spikes up to time t. It follows from (1) that the probability of a single spike in a small interval (t, t + Δ] is approximately
equation M2
(2)
Details can be found in [15] and [34]. When A. is small, (2) is approximately the spiking propensity at time t.
The CIF generalizes the rate function of a Poisson process to a rate function that is history dependent. Because the conditional intensity function completely characterizes a spike train, defining a model for the CIF defines a model for the spike train [12], [13]. For our analyses, we use the GLM to define our CIF models by expressing for each neuron, the log of its CIF in terms of the neurons spike history and relevant movement covariates [35]. The GLM is an extension of the multiple linear regression model in which the variable being predicted, in this case spike times, need not be Gaussian [27]. GLM also provides an efficient computational scheme for model parameter estimation and a likelihood framework for conducting statistical inferences [12].
We express the CIF for each neuron as a function of movement direction, which corresponds to up, right, left, and down, and the neuron’s spiking history in the preceding 150 ms. Instead of estimating the CIF continuously throughout the entire trial, we estimate it over 350 ms time windows around key epochs and at discrete time intervals each 1 ms in duration. Specifically, we estimate the CIF over 350 ms windows centered at the gray array onset (GA), TC onset, GC onset, and movement onset (MV) onsets. Fig. 2 shows all of the time periods for which we estimate the CIF. Henceforth, we omit the superscripts denoting the epoch for a simpler read and express the CIF as follows:
equation M3
(3)
where λS (t|Θ) describes the effect of the movement direction stimulus on the neural response and λH (t|Ht, Θ) describes the effect of spiking history on the neural response. Θ is a parameter vector to be estimated from data. The units of λS(t|Θ) is spikes per second and λH (t|Ht, Θ) is dimensionless. The idea to express the CIF as a product of a stimulus component and a temporal or spike history component was first suggested in [25] and is appealing as it allows one to assess how much each component contributes to the spiking propensity of the neuron. If spiking history is not a factor associated with neural response, then λH (t|Ht, Θ) will be very close to 1 for all times and (3) reduces to an inhomogeneous Poisson process.
Fig. 2
Fig. 2
Time periods, over which the CIF denoted by (3) is estimated, are shaded.
The model of the stimulus effect is as follows:
equation M4
(4)
where
equation M5
The equation M6 parameters measure the effects of movement direction on the spiking propensity. For example, if α1 is significantly larger than α2, α3, and α4 during movement, then the probability that a neuron will spike is greater when the patient moves in the up direction, suggesting that the neuron may be tuned in the up direction.
Our model of spike history effect is as follows:
equation M7
(5)
where n(a : b) is the number of spikes observed in the time interval [α, b) during the epoch. The equation M8 parameters measure the effects of spiking history in the previous 10 ms, and therefore, can capture refractoriness and/or bursting on the spiking probability in the given epoch. For example, if eβ1 is close to zero for any given epoch, then for any given time t, if the neuron had a spike in the previous millisecond then the probability that it will spike again is also close to zero (due to refractory period). Or if eβ5 is significantly larger than 1, then for any time t, if the neuron had a spike five prior to t, then the probability that it will spike again is modulated up, suggesting bursting.
The equation M9 parameters measure the effects of the spiking history in the previous 10–150 ms on the spiking probability, which may be associated with not only the neuron’s individual spiking activity, but also that of its local neural network. For example, if eγ4 is significantly larger than 1, then for any time t if the neuron had one or more spikes between 40–50 prior to t, then the probability that it will spike again is modulated up, suggesting 20–25 Hz oscillations.
By combining (4) and (5), we see that the CIF may be written as follows:
equation M10
(6)
The model parameter vector Θ = αb, βj, γk contains 28 unknown parameters (for each epoch and time window modeled). We compute maximum-likelihood (ML) estimates for Θ and 95% confidence intervals of Θ for each neuron using glmfit.m in MATLAB.
Finally, to choose the model defined by (6), we varied the history bins [a and b in n(a:b)] and also varied the basis functions Id (t) multiplying the parameters {αd} and analyzed the tradeoff between the number of parameters in the model and the ML cost of the fitted model. That is, we computed the Aikaike’s criterion (AIC) [3], which is 2*(number of parameters – log likelihood of model), for each possible model. The model given by (6) was selected as optimal because it rendered the minimum AIC amongst model classes explored.
F. Model Fitting
Establishing the degree of agreement between a point process model and observations of the spike train and associated experimental variables is a prerequisite for using the point process analysis to make scientific inferences. We used Kolmogorov–Smirov (KS) plots based on the time-rescaling theorem to assess model goodness-of-fit. The time-rescaling theorem is a well known result in probability theory, which states that any point process with an integrable conditional intensity function may be transformed into a Poisson process with unit rate [24]. A KS plot, which plots the empirical cumulative distribution function of the transformed spike times versus the cumulative distribution function of a unit rate exponential, is used to visualize the goodness-of-fit for each model. The model is better if its corresponding KS plot lies near the 45° line. We computed 95% confidence bounds for the degree of agreement using the distribution of the KS statistic [24]. If a model’s KS plot was within the 95% confidence bounds, we included it in our analyses.
As mentioned earlier, we built point process models for STN neurons in seven Parkinson’s patients and one healthy primate, which captured dynamics across four different epochs within a directed hand-movement task. We summarize results for each species later. For the PD data, 28 STN neuron models passed the KS test and for the primate data, 35 models passed the KS test.
Recall from (2) that λ(t|Ht)Δ is approximately the probability that the neuron will spike at time t given extrinsic and intrinsic dynamics up to time t, which is captured in Ht. By virtue of (6), we allow the probability that each STN neuron will spike at some time t to be modulated by movement direction, short-term history, and long-term history spiking dynamics. Fig. 3 illustrates these three modulation factors on spiking activity for both PD and primate single neuron models by plotting the optimal parameters and their corresponding 95% confidence bounds before and after MV onset. We make the following observations.
Fig. 3
Fig. 3
Optimal model parameters for an STN neuron during MV– and MV+ periods of a (left) PD patient and (right) healthy primate. [Top row (movement direction modulation)] Optimal extrinsic factors eαd for d = 1,2,3,4(11, R, D, L) are plotted (more ...)
  • Refractoriness: As illustrated in the second row of Fig. 3, both the PD and primate STN neuron exhibits refractory periods [9] as indicated by down modulation by a factor of ten or more due to a spike occurring 1 ms prior to a given time t. That is, if a spike occurs 1 ms prior to time t, then it is very unlikely that another spike will occur at time t (eβ1 < 1 for all eβ1 within its 95% confidence band).
  • Bursting: As illustrated in the second row of Fig. 3, the PD STN neuron spikes in rapid succession before and after MV onset as indicated by one or more of the short-term history parameters (eβi ‘s) corresponding to 2–10 ms in the past being larger than 1. That is, if a spike occurs 2–10 ms prior to time t, then it is more likely that another spike will occur at time t. Formally, a neuron bursts if its model parameters satisfy the following: for at least one i = 2,3,…, 10, LBi ≥ 1 and UBi ≥ 1.5, where LBieβi ≤ UBi, LB and UB are the 95% lower and upper confidence bounds, respectively.
  • 10–30 Hz oscillations: As illustrated in the third row of Fig. 3, the PD STN neuron exhibits 10–30 Hz oscillatory firing before movement. That is, the probability that the PD STN neuron will spike at a given time t is modulated up if a spike occurs 30–100 ms prior to t. Formally, a neuron has 10–30 Hz oscillations if its model parameters satisfy the following: for at least one i = 2, 3,…, 5, LBi ≥ 1 and UBi ≥ 1.5, where LBieγi ≤ UBi.
  • Directional tuning: As illustrated in the first row of Fig. 3, the PD STN neuron appears to exhibit more directional tuning after MV onset. That is, the PD neuron seems more likely to spike in one direction more than at least one other direction. To quantify directional tuning, we performed the following test for each neuron, each time relative to onset, and each epoch:
    • For each direction d* = {U, R, D, L}, compute Pd*d = Pr(eαd* >eαd) = Pr(αd. > αd) for d ≠ d*. Define pd*d* = 0. Use the Gaussian approximation for αd, which is one of the asymptotic properties of ML estimates to compute pd*d.
    • If maxd* = 1,2,3,4 Pd*d ≥ 0.975 then neuron exhibits directional tuning.
Table II provides a population summary for each of these spiking characteristics for each epoch and subject group.
TABLE II
TABLE II
Population Summary for PD Patients and Primate
Fig. 4 plots the population summary for each subject group and also marks with a “+” sign when the population summary for a given characteristic is statistically significantly above or below (with a p-value ≤ 0.05) the population summary baseline value, which is taken as right before the GA appears (GA–). Sign tests are used to test the null hypothesis that the median of a distribution (in our case, the median is baseline) is equal to some value [39]. Under the null hypothesis, we expect half of the observations to be above the median and half to be below, therefore, the number of observations at any given time window during a trial (e.g., MV– is a 350 ms window right before MV onset) that are above the median, which follow a binomial distribution with p = 1/2 and n equals the total number of observations above the mean and below the mean. We considered each patient separately and computed the number of observations (across all neurons in patient) that collectively were above or below the patient’s baseline for each STN characteristic (bursts, high-frequency oscillations (HFOs), and directional tuning) and for each epoch. We then computed either 1, if the probability of observing something greater if the number of observations were greater than the baseline under the null hypothesis, or 2, if the probability of observing something less if the number of observations were less than the baseline under the null hypothesis. These probabilities are the p-values.
Fig. 4
Fig. 4
STN population summary using point process model parameters. Dashed line: bursting, dotted line: 10–30 Hz oscillations, and solid line: directional tuning.
We make the following observations from Fig. 4 and Table II. Most neurons in both species exhibit refractoriness. Bursting is prevalent across all epochs in neural activity of PD patients (on average 39% of PD STN neurons burst). In contrast, neural activity in the healthy primate exhibits little bursting (14% on average) across all epochs. Ten to thirty hertz oscillations are prevalent in neural activity of PD patients during across all epochs (on average 36%) and significantly decrease relative to this baseline post movement as denoted by “+” symbols in the solid curve at the top of Fig. 4. Beta oscillations have been observed experimentally in both Parkinsonian primates and PD patients [7], [8], [10], [17], [29], [32], and attenuation of these oscillation post movement have also been observed [4], [36]. In contrast, an average of 12% of the primate neurons exhibit 10–30 Hz oscillations, which does not significantly modulate across the entire trial. Directional tuning is more prevalent in the healthy primate across the trial. In particular, directional tuning increases significantly above baseline right after the GA is shown in the primate case (see “+” in solid curve at the bottom of Fig. 4 at GA+). This makes sense as the primate knows and moves to one of the four possible directions shown. Tuning increased further in the primate neurons after the TC appears, as now the subject knows which direction to move when cued to move. In contrast, directional tuning fails to increase significantly above baseline until right before MV onset (see “+” in solid curve at the top plot of Fig. 4 at MV–) in PD STN neurons. The lack of significant increase in directional tuning in PD STN neurons early on in the trial may reflect the lack of a dynamic range in the STN neurons of PD patients, which may cause their slow and impaired movements.
For comparison, we also computed spiking characteristics using traditional methods. Next, we describe these computations.
A. Beta Oscillations
To analyze beta oscillations, we computed spectrograms for each epoch window (e.g., 350 ms window right before the GA appears), each neuron and each trial. Then, for each spectrogram, we computed the oscillation density of the spectrogram in the 10–30 Hz range as the integral of the spectrogram in the 10–30 Hz range divided by the integral of the spectrogram across all frequencies. That is, both were double integrals computed across specified frequencies in the 10–30 Hz range and across all time samples in the epoch. Then, for a given neuron and a given epoch, we computed the fifth percentile across all trials as a lower confidence bound on oscillation density (LBOS). We determined the neuron as exhibiting 10–30 Hz oscillations if LBOS > 0.155.
B. Bursting
To analyze bursting in these neurons, we computed interspike interval (ISI) histograms across each epoch during the trial for all neurons and all trials. We then normalized each histogram, so that it is summed to 1, and then, computed the bursting density of the histogram in the 2–10 ms range by taking the sum of the normalized histogram in the 2–10 ms range. Once we computed densities across all trials for a given neuron and epoch, we computed the fifth percentile as a lower confidence bound (LBBU). For a given neuron and epoch, we determined that the neuron bursts if the LBBU > 0.15.
C. Directional Tuning
To analyze directional tuning in these neurons, we computed tuning vectors [20], [21], [33] across each epoch during the trial for all neurons and all trials. If the vector sum in all four directions lies within 20° from one of the four directions, we determined that the neuron is directionally tuned.
The population summary using traditional statistics is shown in Fig. 5 for each subject group. When comparing Figs. 4 and and5,5, we see similar trends in the spiking characteristics of primate STN neurons though absolute percentages slightly differ. In particular, we see that for the primate, we have a steady average of 17% neurons bursting and 19% neurons oscillating in the 10–30 Hz range over the entire trial. We also see significantly increased directional tuning relative to baseline right after the GA appears (GA+). One visible difference between the two analyses (in the primate case) is that the point process models show directional tuning continuously increasing after the TC is shown (solid curve at the bottom of Fig. 4), whereas, the tuning vector analysis shows a decrease in directional tuning around the GC epoch (solid curve at the bottom of Fig. 5). The reason for this discrepancy is due to the fact that tuning vectors are only capturing first-order statistics of the point process. The point process model parameters take into account the stimulus parameters (α’s) probability distributions (not just the mean values), and directional tuning is determined from these distributions as described earlier. Therefore, making inferences from average tuning vectors can be misleading.
Fig. 5
Fig. 5
STN population summary using traditional statistics. Dashed line: bursting, dotted line: 10–30 Hz oscillations, and solid line: directional tuning.
In the human case, we see significant differences between the two analyses. The point process models show significant bursting and oscillations throughout the trial (an average of 39% and 36%, respectively), while traditional methods lead us to believe that there is much less bursting than 10–30 Hz oscillations (an average of 10% and 37%, respectively). The reason the ISI histogram does not show as much bursting in the neurons is precisely because there are also 10–30 Hz oscillations in the spiking activity. Therefore, there are secondary peaks in the ISI histograms between 30–100 ms range. These secondary peaks result in less bursting density, leading us to believe that bursting may not be prevalent. The fact is, PD STN neurons often burst and oscillate in the beta frequency range and this is captured by the point process model parameters as described earlier, since our GLM separates the contributions of short-term intrinsic factors and long-term intrinsic factors on the propensity of the neuron to spike (5).
Another drawback of using traditional statistics is that they are significantly different from one another, and therefore, using them to define whether a neuron bursts, oscillates, or exhibits directional tuning over a certain epoch is not straightforward. In fact, the population summary shown in Fig. 5 varies significantly as we change thresholds. In contrast, for point process models, we can use the same threshold to determine whether a neuron oscillates and bursts as the threshold is on how model parameters modulate the overall probability that the neuron spikes at any given time.
We have applied the point process framework to the analysis of STN microelectrode recordings from PD patients and a healthy nonhuman primate, to understand the relative importance of movement and spiking history on neural responses. We used GLM representations of the point process CIF to develop an efficient likelihood-based approach to model fitting, goodness-of-fit assessment, and inference. The point process model parameters allowed us to identify pathological characteristics of the STN neurons in PD patients, including bursting, 10–30 Hz oscillations, and decreased directional tuning prior to movement. These characteristics, which differed from the characteristics of the non-PD STN neurons, had been previously described using traditional methods. However, such techniques can lead to erroneous inferences when spiking data contains significant temporal dependencies as is the case of PD STN spiking activity. The point process framework is, therefore, a useful paradigm for providing a succinct, quantitative characterization of the pathological behavior of STN spiking activity in PD patients.
Acknowledgments
The work of S. V. Sarma was supported by Burrough’s Wellcome Fund and a L’Oreal FWIS Fellowship. The work of E. Eskandar was supported in part by NIH—National Eye Institute under Grant 1R01EY017658, in part by NIH—National Institute Drug Addiction under Grant 1R01NS063249, in part by National Science Foundation under Grant IOB 0645886, and in part by HHMI and Klingenstein Foundation. The work of E. N. Brown was supported under Grants R01 DA015644 and DP1 OD003646.
Biographies
An external file that holds a picture, illustration, etc.
Object name is nihms506198b1.gif Object name is nihms506198b1.gif
Sridevi V. Sarma (M’04) received the B.S. degree in electrical engineering from Cornell University, Ithaca, NY, in 1994, and the M.S. and Ph.D. degrees in electrical engineering and computer science from Massachusetts Institute of Technology, Cambridge MA, in 1997 and 2006, respectively.
From 2006 to 2009, she was a Postdoctoral Fellow in the Brain and Cognitive Sciences Department, Massachusetts Institute of Technology. She is currently an Assistant Professor in the Department of Biomedical Engineering, Institute for Computational Medicine, Johns Hopkins University, Baltimore MD. Her research interests include modeling, and estimation and control of neural systems.
Dr. Sarma is the recipient of the General Electric Faculty for the future scholarship, a National Science Foundation Graduate Research Fellow, a L’Oreal For Women in Science Fellow, and the Burroughs Wellcome Fund Careers at the Scientific Interface Award.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b2.gif Object name is nihms506198b2.gif
Uri T. Eden received the B.S. degree in mathematics, and engineering and applied sciences from the California Institute of Technology, Pasadena, CA, in 1999, the S.M. degree in engineering sciences, and the Ph.D. degree in engineering sciences from Harvard University, Cambridge, MA, in 2002 and 2005, respectively.
From 2005 to 2006, he was a Postdoctoral Fellow in the Brain and Cognitive Sciences Department, Massachusetts Institute of Technology, Cambridge. In 2006, he joined the Department of Mathematics and Statistics, Boston University, Boston, MA. His current research interests include developing mathematical and statistical methods to analyze neural spiking activity, integrating methodologies related to model identification, statistical inference, signal processing, and stochastic estimation and control.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b3.gif Object name is nihms506198b3.gif
Ming L. Cheng received the B.S. degree in chemistry, molecular and cellular biology, and nutrition from the University of California, Berkley, CA, in 1994, and the M.D. degree from Stanford Medical School, Palo Alto, CA, in 1998.
From 1999 to 2001, he was a Resident in Neurosurgery, Barrow Neurological Institute, Phoenix, AZ. From 2002 to 2005, he was in Neurosurgery, University of North Carolina at Chapel Hill, Chapel Hill, NC. From 2006 to 2008, he was a Clinical and Research Fellow in Functional and Epilepsy Neurosurgery, Harvard Medical School, Massachusetts General Hospital, Boston, MA. He is currently an Assistant Professor in the Department of Neurosurgery, Alpert Medical School, Brown University, Providence, RI, and a Chief of Functional and Epilepsy Surgery, Rhode Island Hospital/The Neurosurgery Foundation, Providence. His research interests include deep brain stimulation and epilepsy.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b4.gif Object name is nihms506198b4.gif
Ziv M. Williams received the B.S. degree in biochemistry and cell biology from the University of California, San Diego, CA, in 1994, and the M.D. degree from Stanford Medical School, Palo Alto, CA. in 1999.
From 1998 to 1999, he was a Doctoral Research Fellow in the Department of Neurosurgery, Stanford School of Medicine. From 1999 to 2000, he was a General Surgery Intern at Harvard Medical School, Massachusetts General Hospital, Boston, MA, where he is currently a Postdoctoral Research Fellow in the Department of Neurobiology. From 2000 to 2006, he was a Resident in Neurosurgery and from 2003 to 2004, a Fellow in Stereotactic and Functional Neurosurgery at Massachusetts General Hospital.
Dr. Williams is a member of the Society for Neuroscience, the American Association of Neurological Surgeons, and the American Society of Stereotactic and Functional Neurosurgery.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b5.gif Object name is nihms506198b5.gif
Rollin Hu received the B.S. and A.B. degrees in biological sciences and public policy from Stanford University, Stanford, CA, in 2000, and the M.D. degree from Harvard Medical School, Boston, MA, in 2004.
In 2005, he was a General Surgery Intern, where he is currently a Resident in Neurosurgery, Massachusetts General Hospital (MGH), Boston. He has been a Postdoctoral Research Fellow at the Department of Neurosurgery, MGH-Harvard Medical School Center for Nervous System Repair.
Dr. Hu is a member of the American Association of Neurological Surgeons and the Congress of Neurological Surgeons.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b6.gif Object name is nihms506198b6.gif
Emad Eskandar received the B.A. degree in chemistry from the University of Nebraska, Lincoln, NE, in 1987, and the M.D. degree from the University of Southern California, Los Angeles, CA, in 1993.
From 1989 to 1991, he was a Fellow in Neurophysiology, National Institutes of Health, Bethesda, MD. From 1993 to 1999, he was a Resident in Neurosurgery, Massachusetts General Hospital, Boston, MA. From 1996 to 1998, he was a Postdoctoral Fellow in Neurophysiology, Harvard Medical School, Cambridge, MA. He is currently a Neurosurgeon at the Massachusetts General Hospital and an Associate Professor of Surgery at Harvard Medical School. He has been involved in multidisciplinary program and also in psychiatry, neurology, and anesthesia. He is also engaged in cuttingedge surgical treatments for movement disorders, epilepsy, pain, and severe psychiatric disorders. He has been involved in basic research laboratory, where he is engaged in theories of learning, motivation, depression and drug addiction.
An external file that holds a picture, illustration, etc.
Object name is nihms506198b7.gif Object name is nihms506198b7.gif
Emery N. Brown (M’01–SM’06–F’08) received the B.A. degree from Harvard College, Cambridge, MA, the M.D. degree from Harvard Medical School, Boston, MA, and the A.M. and Ph.D. degrees in statistics from Harvard University, Cambridge, MA.
He is currently a Professor of computational neuroscience and health sciences and technology in the Department of Brain and Cognitive Sciences and the Division of Health Sciences and Technology, Harvard-Massachusetts Institute of Technology, the Warren M. Zapol Professor of Anesthesia at Harvard Medical School, and an Anesthesiologist in the Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital. His research interests include signal processing algorithms for the study of neural systems, functional neuroimaging, and electrophysiological methods to study general anesthesia in human subjects.
Dr. Brown is a member of the Association of University of Anesthesiologists. He is also a Fellow of the American Statistical Association, a Fellow of the American Association for the Advancement of Science, and a member of the Institute of Medicine. He was a recipient of the 2007 National Institutes of Health Director’s Pioneer Award.
Contributor Information
Sridevi V. Sarma, Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
Uri T. Eden, Mathematics and Statistics Department, Boston University, Boston, MA 02115 USA.
Ming L. Cheng, Clinical Neuroscience Department, Brown University, Providence, RI 02912 USA.
Ziv M. Williams, Neurosurgery, Massachusetts General Hospital, Boston, MA 02114 USA.
Rollin Hu, Neurosurgery, Massachusetts General Hospital, Boston, MA 02114 USA.
Emad Eskandar, Neurosurgery, Massachusetts General Hospital, Boston, MA 02114 USA.
Emery N. Brown, Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
1. Abosch A, Hutchison WD, Saint-Cyr JA, Dostrovsky JO, Lozano AM. Movement-related neurons of the subthalamic nucleus in patients with Parkinson disease. J. Neurosurg. 2002;vol. 97:1167–1172. [PubMed]
2. Albin RL, Young AB, Penney JB. The functional anatomy of basal ganglia disorders. Trends Neurosci. 1989;vol. 12:366–375. [PubMed]
3. Akaike H. A new look at the statistical model identification. IEEE Trans. Autom. Control. 1974 Dec;vol. AC-19(no. 6):716–723.
4. Arminovin R, Williams ZM, Cosgrove GR, Eskandar EN. Visually guided movements suppress subthalamic oscillations in Parkinson’s disease patients. J. Neurosci. 2004;vol. 24(no. 50):11302–11306. [PubMed]
5. Arminovin R, Williams ZM, Cosgrove GR, Eskandar EN. Experience with microelectrode guided subthalamic nucleus deep brain stimulation. Neurosurgery. 2006 Feb;vol. 58(no. 1) [PubMed]
6. Barbieri R, Quirk MC, Frank LM, Wilson MA, Brown EN. Construction and analysis of non-poisson stimulus response models of neural spike train activity. J. Neurosci. Methods. 2001;vol. 105:25–37. [PubMed]
7. Bergman H, Wichman T, Karmon B, DeLong MR. The primate subthalamic nucleus. II. Neuronal activity in the MPTP model of parkinsonism. J. Neurophysiol. 1994;vol. 72:507–520. [PubMed]
8. Bevan MD, Magill PJ, Terman D, Bolam IP, Wilson CI. Move to the rhythm: Oscillations in the subthalamic nucleus-external globus pallidus network. Trends Neurosci. 2002;vol. 25:525–531. [PubMed]
9. Brodal P. The Central Nervous System: Structure and Function. New York: Oxford Univ. Press; 1998.
10. Brown P. Oscillatory nature of human basal ganglia activity: Relationship to the pathophysiology of Parkinson’s disease. Mov. Disord. 2003;vol. 18:357–363. [PubMed]
11. Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM. The time-rescaling theorem and its application to neural spike train data analysis. Neural Comput. 2002;vol. 14:325–346. [PubMed]
12. Brown EN, Barbieri R, Eden UT, Frank LM. Likelihood methods for neural data analysis. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach. London: CRC; 2003. pp. 253–286. ch. 9.
13. Brown EN. Theory of point processes for neural systems. In: Chow CC, Gutkin B, Hansel D, Meunier C, Dalibard J, editors. Methods and Models in Neurophysics. ch. 14. Paris: Elsevier; 2005. pp. 691–726.
14. Brown EN, Kass RE, Mitra P. Multiple neural spike train data analysis: State-of-the-art and future challenges. Nature Neurosci. 2004;vol. 7:456–161. [PubMed]
15. Cox DR, Isham V. Point Processes. Boca Raton, FL: CRC; 2000.
16. Crutcher MD, DeLong MR. Single cell studies of the primate putamen. II. Relations to direction of movement and pattern of muscular activity. Exp. Brain Res. 1984;vol. 53(no. 2):244–258. [PubMed]
17. Dostrovsky I, Bergman H. Oscillatory activity in the basal ganglia— Relationship to normal physiology and pathophysiology. Brain. 2004;vol. 127:721–722. [PubMed]
18. Eden UT, Amirnovin R, Brown EN, Eskandar EN. Proc. Joint Stat. Meet. (JSM) Salt Lake City, UT: 2007. Constructing models of the spiking activ of neurons in the subthalamic nucleus of Parkinson’s patients.
19. Gdowski MJ, Miller LE, Parrish T, Nenonene EK, Houk JC. Context dependency in the globus pallidus internal segment during targeted arm movements. J. Neurophysiol. 2001 Feb;vol. 85(no. 2):998–1004. [PubMed]
20. Georgopoulos AP, DeLong MR, Crutcher MD. Relations between parameters of step-tracking movements and single cell discharge in the gobus pallidus and subthalamic nucleus of the behaving monkey. J. Neurosci. 1983;vol. 3:1586–1598. [PubMed]
21. Georgopoulos AP, Kalaska JF, Caminiti R, Massey JT. On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. J. Neurosci. 1982;vol. 2:1527–1537. [PubMed]
22. Hudsen D, Cohen ME. Chapter in Neural Engineering: Neural Signal Processing. New York: Springer-Verlag; 2005. pp. 193–219.
23. Hutchison WD, Allan RJ, Opitz H, Levy R, Dostrovsky JO, Lang AE, Lozano AM. Neurophysiological identification of the subthalamic nucleus in surgery for Parkinson’s disease. Ann. Neurol. 1998;vol. 44(no. 4):622–628. [PubMed]
24. Johnson NL, Kotz S. Continuous Univariate Distributions–2. New York: Wiley; 1970. Distributions in statistics.
25. Kass RE, Ventura V. A spike train probability model. Neural Comput. 2001;vol. 13:1713–1720. [PubMed]
26. Levy R, et al. Dependence of subthalamic nucleus oscillations on movement and dopamine in Parkinson’s disease. Brain. 2002;vol. 125:1196–1209. [PubMed]
27. McCullagh P, Nelder JA. Generalized Linear Models. 2nd ed. Boca Raton, FL: Chapman and Hall/CRC; 1989.
28. Mitra PP, Jarvis MR. Sampling properties of the spectrum and coherency of sequences of action potentials. Neural Comput. 2001 Apr;vol. 13(no. 4):717–749. DOI: arXiv:physics/0002053vl. [PubMed]
29. Montgomery E., Jr Subthalamic nucleus neuronal activity in Parkinson’s disease and epilepsy subjects. Parkinsonism Relat. Disord. 2008;vol. 14(no. 2):120–125. [PubMed]
30. Paninski L. Maximum likelihood estimation of cascade point-process neural encoding models. Netw.: Comput. Neural Syst. 2004;vol. 15(no. 4):243–262. [PubMed]
31. Paradiso G, Saint-Cyr JA, Lozano AM, Lang AE, Chen R. Involvement of the human subthalamic nucleus in movement preparation. Neurology. 2003;vol. 61(no. 11):1538–1545. [PubMed]
32. Raz A, Vaadia E, Bergman H. Firing patterns and correlations of spontaneous discharge of pallidal neurons in the normal and the tremulous l-methyl-4-phenyl-l,2,3,6-tetrahydropyridine vervet model of parkinsonism. J. Neurosci. 2000;vol. 20:8559–8571. [PubMed]
33. Schwartz AB. Direct cortical representation of drawing. Science. 1994;vol. 265:540–542. [PubMed]
34. Snyder DL, Miller MI. Random Point Processes in Time and Space. New York, NY: Springer-Verlag; 1991.
35. Truccollo W, Eden UT, Fellow MR, Donoghue JP, Brown EN. A point process framework for relating neural spiking activity for spiking history, neural ensemble and extrinsic covariate effects. J. Neurophys. 2005;vol. 93:1074–1089. [PubMed]
36. Williams ZM, Neimat JS, Cosgrove GR, Eskandar EN. Timing and direction selectivity of subthalamic and pallidal neurons in patients with Parkinson disease. Exp. Brain Res. 2005 May;vol. 162(no. 4):407–116. [PubMed]
37. Brody CD. Correlations without synchrony. Neural Comput. 1999;vol. 11:1537–1551. [PubMed]
38. Keat J, Reinagel P, Reid RC, Meister M. Predicting every spike: A model for the responses of visual neurons. Neuron. 2001;vol. 30(no. 3):803–817. [PubMed]
39. Zar JH. Biostatistical Analysis. 4th ed. Upper Saddle River, NJ: Prentice-Hall; 1999.
40. Zelnikera E, Bradley AP, Castnerc JE, Cheneryc HJ, Copland DA, Silburnd PA. Estimation of neuronal firing rates with the three-state biological point process model. J. Neurosci. Methods. 2008 Sep 30;vol. 174(no. 2):281–291. [PubMed]
41. Hoehn M, Yahr M. Parkinsonism: Onset, progression and mortality. Neurology. 1967;vol. 17(no. 5):427–142. [PubMed]
42. Czanner C, Eden UT, Wirth S, Yanike M, Suzuki W, Brown EN. Analysis of between-trial and within-trial neural spiking dynamics. J. Neurophysiol. 2008;vol. 99:2672–2693. [PMC free article] [PubMed]
43. Truccolo W, Hochberg LR, Donoghue JP. Collective dynamics in human and monkey sensorimotor cortex: Predicting single neuron spikes. Nat. Rev. Neurosci. 2010;vol. 11(no. 2) Nat. Neurosci. (2009, Dec.6).[Online]. Doi: 10.1038/nn.2455. (News on Research Highlights) [PMC free article] [PubMed]