--------------------
GENERAL INFORMATION 
--------------------
This readme file was generated on 2024-08-27 by Daisuke Kase in Turner lab at the University of Pittsburgh.

Title of Dataset: MatLab Code for Analyzing Reach-related Single Unit Activity in the Internal Segment of the Globus Pallidus in Parkinsonian Macaque
Description of Dataset: This dataset includes MatLab codes for analyzing the electrophysiological data collected in Turner lab at the University of Pittsburgh.
Principal Investigator: Robert S. Turner, rturner@pitt.edu, https://orcid.org/0000-0002-6074-4365
Date of Data Collection: 2015-03-03 through 2019-01-31
Software Dependencies: All codes are written in MatLab (MathWorks, RRID: SCR_001622). In addition to MatLab, the MatNWB interface (Neurodata Without Borders, RRID: SCR_021156) is required to convert our NWB format data on the DANDI Archive (DANDI:000947/0.240510.2211) to the MatLab format data.
Update History: 2025-04-24 Two changes were made by Daisuke Kase
                 1. "make_figures_from_source_data.m" in the "Code/Step4_Figure_and_Table" directory was modified to plot the figures in the revised manuscript.
                 2. "make_table_from_source_data.m" in the "Code/Step4_Figure_and_Table" directory was debugged to generate a correct table.


-------------
FILE OVERVIEW 
-------------
Directory of Files: 
GPi_MPTP/
	_README_general.txt
	Code/
		_README_code.txt
		Step1_Reconstruct_Database/
				NWB2Mat.m
		Step2_Process_Data/
				BasicSpikeMetrics.m
				OscillationAnaly.m
				PeriEventSDF.m
				PeriMvtTimeLock.m
				PeriRespPeakSDF.m
				src/
		Step3_Summarize_in_TabularData/
				make_source_data_mvt.m
				make_source_data_rest.m
				src/
		Step4_Figure_and_Table/
				make_figures_from_source_data.m
				make_table_from_source_data.m
				src/

Relationship Between Files: The "Code" directory contains all our MatLab codes stored separately in four sub-directories for each step of the analysis. 

-----------
METHODOLOGY
-----------
Description of methods used for data collection: 
note: This methodology section is continued from the methodology section in the _README_general.txt file.

(5) Reconstruction of data set for the analyses on MatLab (relevant code: NWB2Mat.m)
   Since our analysis codes are written in MatLab, we need to download the entire dataset from the DANDI Archive (DANDI:000947/0.240510.2211) and convert the NWB format data into MATLAB format using the NWB2Mat code. This code will generate a series of datasets for each individual recording session, along with lists of GPi neurons required to run our analysis scripts.

(6) Data processing and analysis
 6-a. Discharge rate and pattern at rest (relevant code: BasicSpikeMetrics.m and OscillationAnaly.m)
   The effects of MPTP intoxication on resting neuronal activity were quantified using standard methods (1, 2). Those analyses were applied to spike trains during the start-position hold periods (SPHPs) of the behavioral task  intervals of unpredictable duration (2-10 sec) during which an animal held its left hand stationary at the start position while awaiting onset of the go-cue. A neurons mean firing rate was calculated as the total number of spikes detected across all hold periods divided by the summed duration of all hold periods. Episodes of burst firing  discrete period of markedly elevated firing rate  were detected using the Poisson Surprise Method (2, 3). Bursts were defined as groups of 4 or more spikes whose inter-spike intervals (ISIs) were unusually short compared with other ISIs of the same spike-train. We used a surprise threshold of 5, which equates to p<0.05 that the candidate burst would occur as a part of a Poisson-distributed sequence of spikes. The overall burstiness of a cell was quantified as the fraction of spikes that occurred during bursts relative to the total number of spikes found across all SPHPs. 
   Rhythmic modulations in firing rate were detected using a shuffled normalization method (4-6). The discrete Fourier transform (FFT) was applied to non-overlapping 512-ms long segments of a spike trains delta function smoothed with a Hanning window of the same length. (These spike train segments were extracted separately from each SPHP to avoid possible anomalies due to data discontinuities at SPHP concatenation boundaries. SPHPs containing fewer than four spikes were excluded from this analysis.) The resulting primary spectral density estimate (0500-Hz, 2-Hz resolution) was normalized by dividing it by a control spectrum. The control spectrum was the mean of 100 spectra computed after 100 shufflings of the same ISIs as used to compute the primary spectrum. This normalization compensated for distortions in spectral estimates attributable to a neurons refractory period and thereby improved the detection of low frequency oscillations (5). The shuffle-normalization procedure yielded spectra that varied around a normalized value of 1. Peaks in the normalized spectrum between 4-Hz and 100-Hz were tested for significance relative to the SD of the spectrum in the 150250-Hz control range. The omnibus threshold for significance (p=0.05) included a Bonferroni correction for multiple comparisons (0.05 / 51 spectral points tested between 4-Hz and 100-Hz; actual threshold p<9.8^10-4). If a spectral peak exceeded the threshold at more than one frequency bin, then the central rhythmic frequency was defined as the spectral bin with the highest power. 

 6-b. Peri-movement discharge (relevant code: periEventSDF.m)
   We tested for peri-movement changes in single-unit spike rate using a method described previously (7). First, spike density functions (SDFs) were constructed by convolving each units spike time stamps (1 kHz resolution) with a Gaussian kernel (sigma = 25 ms). Spike trains were aligned to the time of movement onset in each trial, and across-trial mean SDFs were constructed separately for reaches to each target. Significant modulations in the mean SDF were tested for during a test window that extended from the across-trials median time of go-cue onset to the across-trials median time of movement offset. The significance of modulations during that test window was determined relative to baseline defined as the mean and standard deviation of the mean SDF across a 700 ms window that ended at the beginning of the test window (defined above), after correcting for linear trends in the mean SDF over the baseline window. We defined a peri-event response as a statistically significant elevation or depression in the mean SDF that lasted at least 60 ms relative to baseline (t-test, one-sample versus baseline; omnibus p < 0.001 after Bonferroni correction for multiple comparisons). Significant peri-movement modulations in mean firing rate were classified as increase- or decrease-type responses. As described previously (8), the peri-movement modulation of some single-units was polyphasic, being composed of one or more significant increase and decrease in firing rate during the test window. Those responses were classified according to the sign of the earliest significant modulation  either polyphasic increase/decrease or polyphasic decrease/increase [referred to hereafter as poly (+/-) and poly (-/+), respectively]. 

 6-c. Quantification of response metrics (relevant code: SDFAnaly_dev2.m)
   The onset latency of a response was defined as the time of the earliest significant point in the response (1 ms resolution) relative to the alignment event of the mean SDF. For monophasic increase- and decrease-type responses, the magnitude of a response was defined as the firing rate maximum or minimum, respectively, relative to the units baseline firing rate estimated at the detected time of response onset. The duration of a response was defined as the full-width of the response at half of its maximum change in firing rate [i.e., full-width at half-max (FWHM), or at half-minimum for decrease-type responses]. For polyphasic responses, response magnitude and duration metrics were computed as outlined above, but only for the initial phase of the response.

 6-d. Trial-by-trial detection of response onsets (relevant code: PeriMvtTimeLock.m and PeriRespPeakSDF.m)
   Measures of response magnitude and duration taken from a mean across-trials SDF are susceptible to distortions that depend on how variable the timing of the response is from trial to trial (i.e., its temporal jitter). For example, in a comparison of across-trials mean SDFs, what appears to be reduction in response magnitude and increase in duration could in fact be attributable to a simple increase in the trial-to-trial variability of response timing. In other words, when measurements are taken from a mean across-trials SDF, the net effect of increased trial-to-trial jitter in response timing is indistinguishable from a true reduction in response magnitude and increase in response duration. 
   To resolve this ambiguity, we detected the onsets of neuronal responses on a trial-by-trial basis using a general approach that has been described in detail elsewhere (8, 9). First, only units found to have a significant perimovement response in the mean across-trials SDF were subjected to this analysis. Among those units, the class of response detected in the mean SDF was used to guide the subsequent trial-by-trial analysis. As an example, the approach used for monophasic perimovement increases is described here. Second, perimovement SDFs for individual trials were computed by convolving spike time stamps (1 ms resolution) with a Gaussian kernel (sigma = 25 ms). Third, for each single trial SDF, we prepared a 200-point (200 ms) response kernel consisting of a simple step function in which the initial 100 points were set to the minimum firing rate in the single-trial SDF and the second 100 points are set to the maximum firing rate. Fourth, the quality of fit between the kernel and segments of the single-trial SDF were computed while shifting the kernel at millisecond steps from the time of go-cue onset to the end of movement. The time shift that resulted in the best fit (minimum summed squared error) between kernel and segment of single-trial SDF was taken to be the time of response onset for that single trial. This procedure was repeated for the single-trial SDFs from each valid behavioral trial. For decrease-type responses, we first inverted the single-trial SDF and then processed as described above for increase-type responses. 
   We then used those trial-by-trial measures of response timing to: 1) analyze the temporal linkage of responses to different task events, 2) estimate the influence of parkinsonism on the fundamental variability of response timing, and 3) test for effects of parkinsonism on response metrics (magnitude and duration) after controlling for potentially confounding influences of differences in a responses trial-to-trial temporal jitter. 
   For item 1) above, Neural responses that participate at different stages of the sensory-motor RT process may be distinguished from each other by the variability in their timing relative to the appearance of the visual go-cue and the onset of movement. If a response has a tight temporal linkage to the sensory go-cue, then the time interval between the go-cue and the neuronal response should be constant regardless of the RT (i.e., slope is about 0), whereas the interval between response onset and movement should covary strongly, across trials, with the behavioral RT (with slope is about 1). If, on the other hand, the response is linked tightly to movement initiation, then the interval between the go-cue and the response should covary with RT. Note that discrimination of those two temporal alignments depends critically on the presence of trial-to-trial variability in the behavioral RT (the interval between go-cue appearance and movement initiation). 
   To quantify those timing relationships, we calculated two separate linear regressions each comparing the trial-by-trial covariation of a time interval (go-cuetoresponse or responsetomovement) with the behavioral RT. A significant positive relation (p<0.05) between responsetomove intervals and RTs, in the absence of a significant relation between cuetoresponse intervals and RTs, was taken to indicate that the response was time locked to go-cue onset (i.e., cue-locked). Conversely, a significant positive relation between cuetoresponse intervals and RTs, but not between responsetomove intervals and RTs, indicated that the response was time locked to movement onset (move-locked). When both regressions yielded significant positive slopes, the response was categorized as having intermediate locking, a phenomenon also reported by (10). If neither regression was significant, then the time locking was categorized as indeterminate. 
   To gain deeper insight into the distribution of different forms of event locking, and the interacting effects of MPTP, we calculated an event locking index (ELI) that reflects the relative difference in slopes yielded by the two regressions (slope of cue-resp and slope of resp-mvt) as follows:  ELI = (slope of cue-resp - slope of resp-mvt)/ (slope of cue-resp + slope of resp-mvt). An ELI value less than 1 denotes perfect cue-locking, ELI more than +1 denotes perfect movement-locking, and an ELI value of 0 denotes a response whose timing covaries with the midpoint of the RT interval regardless of the overall duration of the RT.
   For item 2) above, the fundamental variability of response timing was quantified as the residual trial-to-trial variability in response timing after controlling for potential influences of RT variability and the temporal linkage of the response to specific task events. First, the linear relationship between trial-by-trial response latency and RTs was determined by linear regression (Matlab, FITLM) to model the RT-dependent variability in response timing. Residuals from that regression were taken to reflect the fundamental trial-to-trial variability in response timing. The interquartile range (IQR) of those residuals was used as a measure of the trial-to-trial dispersion in response timing.
   We also used the trial-by-trial estimates of response timing to compare response metrics (magnitude and duration) after controlling for the potentially confounding influences of differences in response temporal jitter. For each response, we removed the effects of temporal jitter by realigning single-trial SDFs to the response itself. First, in individual single-trial SDFs, we found the earliest maximum change in firing rate (Matlab, FINDPEAKS; firing rate maximum for increase-type responses and minimum for decrease-type responses) that occurred following the previously determined time of response onset for that trial. Second, the whole single-trial SDFs were shifted in time so that the times of peak change aligned across trials. Third, those realigned single-trial SDFs were averaged across-trials to yield a response-aligned mean SDF. Finally, response metrics (magnitude and duration) were extracted from the resulting response-aligned mean SDF using methods described above under Quantification of response metrics.

(7) Summarizing the results in tabular format (relevant code: make_source_data_mvt.m and make_source_data_rest.m)
   All results of analyses should be stored in one MatLab format file by running codes described above. Then, make_source_data_rest.m and make_source_data_mvt.m will generate summaries of the results in tabular format. Details of the structure of tabular data are described in _README_source_data.txt file.

(8) Plotting figures and table (relevant code: make_figures_from_source_data.m and make_table_from_source_data.m)
   These codes generate all figures and a table that we use in the next publication. All source data in the SourceData directory are required to run these codes properly.

   

-----------------------
DATA ACCESS AND SHARING
-----------------------
Publications based on this dataset: To be updated
Recommended citation for this dataset: To be updated
License information: To be updated



-----------------------
REFERENCES
-----------------------

1.	Chan V, Starr PA, Turner RS. Bursts and oscillations as independent properties of neural activity in the parkinsonian globus pallidus internus. Neurobiol Dis 41: 2-10, 2011.
2.	Wichmann T, Bergman H, Starr PA, Subramanian T, Watts RL, DeLong MR. Comparison of MPTP-induced changes in spontaneous neuronal discharge in the internal pallidal segment and in the substantia nigra pars reticulata in primates. Experimental brain research Experimentelle Hirnforschung 125: 397-409, 1999.
3.	Legendy CR, Salcman M. Bursts and recurrences of bursts in the spike trains of spontaneously active striate cortex neurons. Journal of neurophysiology 53: 926-939, 1985.
4.	McCairn KW, Turner RS. Deep brain stimulation of the globus pallidus internus in the parkinsonian primate: local entrainment and suppression of low-frequency oscillations. Journal of neurophysiology 101: 1941-1960, 2009.
5.	Rivlin-Etzion M, Ritov Y, Heimer G, Bergman H, Bar-Gad I. Local shuffling of spike trains boosts the accuracy of spike train spectral analysis. Journal of neurophysiology 95: 3245-3256, 2006.
6.	Starr PA, Kang GA, Heath S, Shimamoto S, Turner RS. Pallidal neuronal discharge in Huntington's disease: support for selective loss of striatal cells originating the indirect pathway. Exp Neurol 211: 227-233, 2008.
7.	Zimnik AJ, Nora GJ, Desmurget M, Turner RS. Movement-related discharge in the macaque globus pallidus during high-frequency stimulation of the subthalamic nucleus. J Neurosci 35: 3978-3989, 2015.
8.	Turner RS, Anderson ME. Pallidal discharge related to the kinematics of reaching movements in two dimensions. Journal of neurophysiology 77: 1051-1074, 1997.
9.	Anderson ME, Turner RS. A quantitative analysis of pallidal discharge during targeted reaching movement in the monkey. Experimental brain research Experimentelle Hirnforschung 86: 623-632, 1991.
10.	DiCarlo JJ, Maunsell JH. Using neuronal latency to determine sensory-motor processing pathways in reaction time tasks. Journal of neurophysiology 93: 2974-2986, 2005.

