Delta Oscillations Are a Robust Biomarker of Dopamine Depletion Severity and Motor Dysfunction in Awake Mice

Delta oscillations (0.5–4 Hz) are a robust but often overlooked feature of basal ganglia pathophysiology in Parkinson’s disease and their relationship to parkinsonian akinesia has not been investigated. Here, we establish a novel approach to detect spike oscillations embedded in noise to provide the first study of delta oscillations in awake, dopamine depleted mice. We find that approximately half of neurons in the substantia nigra reticulata exhibit delta oscillations in dopamine depletion and that these oscillations are a strong indicator of dopamine loss and akinesia, outperforming measures such as changes in firing rate, irregularity, bursting and synchrony. We further establish that these oscillations are caused by the loss of D2 receptor activation and do not require motor cortex, contrary to previous findings in anesthetized animals. These results give insight into how dopamine loss leads to dysfunction and suggest a reappraisal of delta oscillations as a biomarker in Parkinson’s disease.

and tend to dissipate under treatments such as dopamine replacement therapy (Weinberger et 23 al., 2006;Ray et al., 2008). Similar oscillations are observed in some animal models of PD -24 slightly higher in frequency (25-35 Hz) in 6-hydroxydopamine (6-OHDA) lesioned rats or lower 25 (8-13 Hz) in 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP)-treated monkeys. In these 26 models, the link between beta oscillations and motor symptoms is less clear. Beta oscillations 27 may arise later than symptoms (Mallet et al., 2008), do not consistently track symptom 28 progression (Muralidharan et al., 2016) or reduction with treatments such as deep brain 29 stimulation (DBS) (McConnell et al., 2012), and occur in both parkinsonian and healthy animals 30 (Connolly et al., 2015). Attempts to artificially induce beta oscillations in these animals have also 31 been insufficient to cause PD-like symptoms (Swan et al., 2019). Even in humans, DBS studies 32 have shown conflicting results between the correlation of beta oscillations and motor symptoms 33 oscillations tend to weaken during stimulation (Kühn et al., 2006(Kühn et al., , 2008 but not in every case 34 (Rossi et al., 2008) or may return before symptoms reemerge (Foffani et al., 2006). 35 In contrast, the lower frequency oscillations observed in human PD patients have been 36 much less studied, despite occurring in a greater number of neurons than beta oscillations in 37 some patients or in the absence of beta oscillations at all (Levy et al., 2002;Du et al., 2018;38 Zhuang et al., 2019). These delta oscillations are often termed "tremor frequency" oscillations 39 due to their typical coherence with Parkinsonian tremor (Bergman et al., 1994), but they may 40 also arise without such coherence (Hurtado et al., 1999). While delta oscillations have an 41 unclear relationship to tremor, their relationship to other PD symptoms such as bradykinesia 42 and rigidity has not been investigated. 43 This lack of attention is surprising, as slower oscillations have also been observed in 44 animal models of PD. In monkeys, oscillations as low as 3-7 Hz have been observed (Raz et  We first sought to reliably quantify these oscillations in dopamine depleted mice. Neural 90 noise is more prevalent in awake than anesthetized animals, and typically manifests in a power 91 law fashion (called "pink" or "flicker" noise) such that it is dominant in low frequencies. Since the 92 oscillations we observe in SNr units in DD were in the range typically tainted by pink noise, we 93 could not reliably detect them using standard approaches based solely on the power spectral 94 density or transformations of it. Specifically, random peaks in the power spectral density atop 95 pink noise, or the pink noise itself, can easily be misidentified as an oscillation of interest (Figure 96 To ensure that delta oscillations were not merely an immune or inflammatory side effect 120 of the injected toxin or cell death, we treated a cohort of animals intraperitoneally with reserpine, 121 a compound that blocks the vesicular monoamine transporter 2 (VMAT2) complex from 122 packaging monoamines into vesicles. This yielded a monoamine (including dopamine) depletion 123 without any intracranial injection or cellular death and produced akinetic symptoms similar to 124 those observed in bilateral 6-OHDA depleted mice. When we recorded three days after the start 125 of daily reserpine injections, these animals exhibited a high proportion of slowly oscillating units 126 in the SNr (33-100% for each animal, 74 of 119 units pooled), similar to bilaterally depleted 127 animals ( Figure 3a). 128 Given the prevalence of beta oscillations in the DD and PD literature, we sought to 129 determine if these animals' SNr units also exhibited beta oscillations. We defined a wide 130 oscillations across humans and common model species (monkey and rat). We saw no increase 132 in the fraction of beta oscillating units after any form of dopamine depletion, with or without our 133 phase shift criterion (Figure 3c- we looked at the relationship between an animal's fraction of units exhibiting a delta oscillation 143 and its level of dopamine neuron loss (as measured by striatal tyrosine-hydroxylase (TH) 144 immunoreactivity). Performing a linear regression to predict %TH remaining from oscillation 145 fraction showed a relatively strong (r 2 = 0.5267) and significant (p < 0.01 from a bootstrapped 146 99% confidence interval, see Methods) relationship between dopamine loss and the fraction of 147 oscillating units (Figure 4a). 148 Since striatal TH immunoreactivity is not a perfect indicator of parkinsonian symptoms, 149 we also used these measures to predict motor behavior. Prior to in vivo recordings, these 150 animals were given a series of behavioral tests to measure their mobility, dexterity, and strength 151   To determine how each parameter informs the model, we shuffled the testing data for 171 that parameter and calculated how much this loss of information increased the MSE of the 172 model (the 'importance" of that parameter). We then estimated 95% confidence intervals for the 173 importance of each parameter (see Methods). The fraction of oscillating units was the only 174 parameter whose confidence interval did not extend below zero (Figure 4e), suggesting that, 175 when the model is built to include oscillations, they are the only parameter that provides reliably 176 predictive information. In other words, while other parameters may provide information, that 177 information is redundant when the fraction of oscillating units is known. 178 To confirm this in another manner, we rebuilt the models using the same cross-validated 179 training and testing sets as above using only a single parameter at a time, or using all of the Using the same procedure as above to predict PC1 of the animals' behavior, we found 188 very similar results to those predicting TH levelsnamely, firing rate, irregularity, burstiness and 189 synchrony provide some information in predicting behavior, particularly when considered The mechanism behind the observed delta oscillations is unclear, but they could arise To quantify this relation, we performed a series of Granger causality regressions, which 249 make no assumptions on the periodicity of the signals. Rather, they simply attempt to predict 250 changes in M1 ECoG based on its own history (the null, autoregressive model) or by 251 additionally including SNr spiking information from a single unit. For each unit, we computed 252 201 separate models predicting M1, each using SNr spiking information at a different lag 253 between -1 (i.e., past spikes) and +1 seconds (i.e., future spikes). Aligning the lag coefficients of 254 the models for a single unit illustrates a periodicity in their values that matches the oscillation 255 in the delta range difficult. Even when it is possible to record from single units, we have 308 demonstrated that low frequency noise can disrupt these spiking signals as well. 309 To reliably detect low frequency spike oscillations in awake animals, we have introduced 310 phase shift as a novel detection technique which utilizes phase information typically discarded 311 from the Fourier transform. Phase shift measures the local stationarity of a signal composed 312 primarily of one frequencya perfect sine wave would have zero phase shift and high power, 313 but a sine wave with a phase that randomly advances would have high phase shift while 314 maintaining high power. This measure can distinguish our signal of interesta single oscillatory 315 signal that shifts in phase only gradually or rarelyfrom low frequency pink noise, a 316 phenomenon that is not restricted to a single frequency, for which phase components measured 317 at individual frequencies may shift rapidly between adjacent windows. 318

Relationship to previous studies on PD oscillations 319
In PD research, much of the oscillation literature has focused on the beta band 320 The low frequency oscillations that we observe resemble those seen in anesthetized 328 mice and rats, although oscillations in awake settings are generally noisier. By performing these 329 experiments in awake mice, this study rules out concerns that oscillations in the basal ganglia

Headbar implantation 399
Animals were anesthetized with 20 mg/kg ketamine and 6mg/kg xylazine and placed in a 400 stereotaxic frame (Kopf Instruments). Anesthesia was maintained throughout surgery with 1.0-401 1.5% isoflurane. All coordinates were measured in mm with AP and ML measured from bregma 402 and DV relative to the dural surface. The scalp was opened and bilateral craniotomies (for later 403 probe insertion) approximately 1.5 x 1.5 mm in size were drilled over SNr (AP: -3.00, ML: 404 ±1.50), GPe (AP 0.00, ML: ±2.12), or STN (AP: -1.70, ML: ±1.52). A custom-made copper or 405 stainless steel headbar was affixed to the mouse's skull with dental cement (Lang Dental). A 406 well of dental cement was then built around the exposed skull (see in vivo recordings) and filled 407 with a silicon elastomer. 408

ECoG connector implantation 424
For experiments involving electrocorticogram (ECoG) recordings, a male gold connector 425 (Ampityco Electronics) was soldered to a stainless steel wire, and the connector was gently 426 lowered above left or right motor cortex (M1, AP: +1.40, ML: ±1.00) such that the wire touched 427 the dural surface then secured in place with dental cement (Lang Dental). 428

Aspiration lesions 429
For experiments involving M1 lesion, a craniotomy was drilled bilaterally over M1 (AP 430 0.0-2.5, ML 1.0:2.5) and the dura was removed. Using a 20-gaugse suction tube (Miltex) 431 attached to a vacuum source, we aspirated cortex to a depth of 2.5 mm across the craniotomy 432 and under portions of the remaining skull, periodically lightly rinsing the area with saline. We 433 filled the lesioned space with triple antibiotic (bacitracin, neomycin, polymyxin) before sealing 434

the craniotomy with a silicon elastomer (Smooth-On) 435
Post-operative care 436 Upon completion of surgery, animals were injected subcutaneously with 0.5 mg/kg 437 ketofen and placed inside their cage half on/half off a heating pad to recover. Dopamine 438 depleted animals were supplied with trail mix and moistened food to maintain weight and 439

Drugs 443
In addition to the drugs used above during surgery, animals were given the following 444 drugs (Sigma-Aldrich, except when specified) dissolved in 0.9% saline (except when specified). to head-fixation for ten minutes, the silicon elastomer was removed and craniotomies were 454 cleaned with saline. Using a micromanipulator (Sutter Instruments), a linear microelectrode 455 probe with sixteen channels spaced 50 µm apart (NeuroNexus) was lowered into the craniotomy 456 at the coordinates listed above for SNr, GPe or STN. After the initial lowering, a ground wire 457 was placed in saline in the dental cement well on the skull. Once the top of the nucleus (SNr, -458 4.0mm, GPe: -3.60mm, STN: -4.00mm from the top of the brain) was found and high firing rate 459 units were observed, the probe was held stable for at least ten minutes prior to recording. waveforms generated a cluster of spikes significantly distinct from other unit or noise clusters (p 490 < .05), 2) the J3-statistic was greater than 1, 3) the Davies-Bouldin statistic was less than 0.5, 491 recording, it was only used in analysis for the time period when its spike cluster satisfied these 493 criteria, and only if its cluster was present for at least three minutes. Data were then imported 494 into MATLAB (MathWorks) in which all further analysis was performed using custom code 495 except when specified. 496 Since units must fire quickly enough to exhibit an oscillation, only units with a firing rate 497 greater than 5 Hz (over 95% of sorted units) were considered for analysis. As ECoG signals 498 were occasionally corrupted for short time windows, generally due to muscle activity, we visually 499 determined a noise threshold for each recording and zeroed any length of signal within 250 500 milliseconds of any data point whose absolute value exceeded that threshold. ECoG signals 501 were then delta (0.5-4 Hz) bandpassed using a 2 nd order Butterworth filter. 502 503

Renewal-Corrected Power Spectrum 505
For each unit, we downsampled its spike train to 1 kHz and split it into segments of 2 12 506 ms, advancing from one segment to the next with time step size = 2 9 ms. For each 507 segment, we calculated its interspike interval (ISI) probability distribution, 0 ( ). We calculated 508

Phase Shift 525
For the kth time segment, we calculated the uncorrected phase ̃ at each frequency: 526 ̃( , ) = tan −1 ( ( ( ))) 527 and made the following correction such that the phase of each frequency is defined relative to 528 the start of the recording rather than the start of the segment: 529 ( , ) = ( + (̃− 2 ( − 1) ), 2 ) − 530 where mod is the modulus operator and is the time step between adjacent segments (here, 531 2 9 /1000 seconds). In other words, for each frequency, imagine a perfect oscillator with zero 532 phase at the start of the recording. For each segment, we determined what phase this oscillator 533 would reach at the start of the segment and defined that phase to be zero for that segment. This 534 correction ensures that a perfect oscillator would have the same corrected phase for every 535 segment. 536 derivative ( , ) by computing the difference of phase across successive time steps and 538 averaged over each difference to obtain the average absolute rate of phase shift: 539 where N is the number of segments. For brevity, we refer to ( ) as the phase shift. 541 542

Oscillation Detection 543
We detected oscillations in a two-step process by first seeking frequencies with high 544 power and then determining whether these frequencies also had low phase shift. 545 To determine whether a unit reached statistically significantly high power at a particular 546 frequency, we found each local maximum of ̂( ), defined as a value higher than its three 547 neighbors on both sides, within the band 0.5-4 Hz (or 7-35 Hz for detecting beta oscillations). 548 We then estimated a 99% confidence interval of renewal-corrected power from the region of 549 ̂( ) between 100 and 500 Hz, correcting for multiple comparisons (Bonferroni correction) of all 550 frequencies in the band of interest. A peak of ̂( ) was considered significant if it fell above this 551 confidence interval. 552 As our second step, we determined if any frequency detected in the previous step had a 553 significantly low phase shift. We estimated a 95% confidence interval of phase shifts from the 554 region of ( ) between 100 and 500 Hz, correcting for multiple comparisons (Bonferroni 555 correction) if multiple frequencies were detected from the PSD. We concluded that an oscillation 556 was present at a frequency with significant power if the phase shift at that frequency fell below 557 this confidence interval. To determine if two units were synchronous, we used the method and parameters 576 outlined in Willard et al. 2019, which determines the fraction of synchronous spikes above 577 chance after correcting for nonstationarity in a unit's firing rate (Willard et al., 2019). In brief, we 578 windowed both spike trains into 12-second segments and zeroed the first and last four seconds 579 of the segment taken from the second spike train. We performed cross-correlation with a 580 maximum lag of four seconds. Since this maximum lag is equal to the length of time zeroed on 581 the second spike train, this ensures a constant number of non-zero-padded comparisons (nc) at 582 each lag, as opposed to traditional cross-correlation in which nc is a function of lag. We divided 583 the cross-correlogram for the segment by the mean value from 0.5-4 seconds on both sides, 584 which allows the correlation's units to be interpreted as the fraction of spikes greater than 585 chance at a given lag (where 1 = chance). We repeated this process on overlapping segments 586 nonstationarity-corrected cross-correlogram. We generated a 99% confidence interval from the 588 data with lag ≥ 0.5 s (which is a reasonable null distribution due to nc, and thus the variance of 589 the correlation estimate, being held constant). We conclude that a pair is synchronous if its 590 normalized cross-correlation at lag zero is larger than the upper boundary of this confidence 591 interval. We built an individual tree on 80% of the data (20 animals) using the fit method of the 618 DecisionTreeRegressor class in the scikit-learn package for Python to predict the percent of TH 619 remaining (Y) from the above neuronal parameters (a set X). In brief, this method places all 620 training data at the topmost node of a tree and calculates the mean squared error (MSE) of this 621 node as if each animal's TH were estimated to be the mean TH of every animal at the node. We 622 determined, for each parameter X, the threshold T that would most reduce the mean squared 623 error (MSE) of the animals if they were to be estimated in two different sets depending on 624 whether their value of X is "greater than" or "less than or equal to" T. We then found the 625 parameter for which the best T most reduces that MSE and split the animals at that node into 626 two new child nodes according to the identified threshold. We iteratively repeated this process 627 at every node until all terminal nodes had two or fewer animals at them, at which point each 628 terminal node is termed a "leaf" of the tree. 629 We tested the remaining 20% of the data (5 animals) using the DecisionTreeRegressor 630 test method, which runs each animal through the tree (picking > or ≤ at each node as 631 determined by the animal's data) until it reaches a leaf. The mean value of Y at each leaf is the 632 prediction for that animal. We computed the error of the tree as the root-mean-squared error 633 (RMSE) of its 5 predictions. 634 We computed a forest of 1000 such trees through subsampling the data into training and 635 testing sets (Monte Carlo cross-validation) and calculated the top and bottom 2.5 percentiles to 636 approximate a 95% confidence interval for the forest. We generated an intercept-only forest 637 (using no parameters in the training set) and oscillation-only forest (using only the fraction of 638 oscillations and an intercept term in the training set) on the same 1000 bootstrapped training 639 and testing sets. 640 The importance of each parameter was determined using a variant on permutation 641 importance. For a given parameter and tree, consider the set S of values for that parameter in 642 the test set. We produced pseudo-test data with every derangement of S (i.e. 5 animals × 44 643 derangements of 5 values = 220 pseudo-test animals with shuffled data for one parameter). The 644 difference between the RMSE of the real test data and the pseudo-test data is the importance of 645 that parameter for that tree. To determine the parameter importance for the entire forest, we 646 approximate a 95% confidence intervals as above from the 1000 trees.

ECoG-Spike Time Series Regression 651
To determine if SNr neurons have a significant lead/lag relationship with M1, we built a 652 series of regression models predicting an M1 ECoG signal from the spiking of a single SNr unit 653 at various lags. First, we binned the ECoG into 10ms bins and defined the dependent variable Y 654 as the difference between adjacent ECoG measurements to reduce nonstationarity. We then 655 built a 10th order autoregressive model of Y which served as the null model. 656 To incorporate SNr firing into the prediction, we calculated the spike density function 657 (SDF) for an SNr unit by convolving its spike train with a Gaussian function with a standard 658 deviation of 100 ms. We then aimed to determine which time shift of the SDF best improves the 659 prediction of the ECoG. One might use a distributed lag model for this task, where the 660 explanatory variables consist of the time shifted ECoG (autoregression) and all considered time 661 shifts of the SNr SDF simultaneously in a single model, but the multicollinearity of the SDF at 662 different time shifts can heavily bias the regression coefficients. Instead we assumed that, if a 663 lag exists by which the unit firing influences the ECoG or vice versa, then there is only one such 664 lag by which this influence occurs. Thus, we could build an individual model for each time shift 665 of the SDF. Each model used the 10 th order autoregressive terms and one SDF term shifted 666 from between -100 and +100 bins (-1000 to +1000 ms) as its explanatory variables. We built 667 201 such models, which covers the entire range of lags at 1 bin increments. 668 To determine if a significant lead/lag existed, we found the best model as determined by 669 its mean squared error (MSE). We then determined if the model at this lag was significantly 670 better than the null autoregressive model by performing an F-test at α < 0.05, correcting for 201 671 comparisons (Bonferroni correction). As choosing ECoG as the independent variable and using 672 autoregressive terms from the past could introduce bias in favor of SNr predicting M1, we also 673 performed these analyses using SNr as the independent variable (i.e. computing a null 674 autoregressive model for SNr spiking and then computing 201 models at distinct ECoG time 675 shifts to compare to the null), and performed the same analysis as above but in backwards time 676 ECoG bandpowers were significantly different across conditions. For comparisons with two 682 groups, a two-sample t-test was performed, unless data were paired before and after a 683 manipulation (e.g. acute drug infusion), in which case a one-sample t-test was performed. For 684 comparisons with multiple groups compared against a control group, a one-way ANOVA was 685 performed, and if this reached significance at the α = 0.05 level, a Dunnett's post-hoc test was 686 performed to determine if there were individual differences comparing groups to control. 687 Asterisks above comparisons in figures correspond to *: p < 0.05, **: p < 0.01, *** p < 0.0001.   overlapping windows (1st row) and its Fourier transform is computed (corrected for the 901 Notice while the peak frequency (red) has consistent phase across windows, an arbitrary noise 905 frequency (blue) has inconsistent phase. We take the absolute circular difference of phases at 906 each frequency (4th row) and compute whether the frequency identified in the power spectrum 907 also has statistically significantly lower phase difference than the control band. A spike train 908 which has both a significant spectral peak and significant phase difference trough at the same 909 frequency is labeled as oscillating. b. Data from two example oscillating units. Top: 910 Autocorrelation exhibiting oscillations. Middle: Significant peaks (red dots) in the PSD 911 surrounded by pink noise. Bottom: The phase difference at these detected frequencies is 912 significantly lower than control frequencies. c. Same as b, but for two units whose 913 autocorrelation appears to be non-oscillating yet have a peak in their PSD and which would be 914 "false positive" detections if only PSD's were analyzed without the consideration of phase shift. Fraction of oscillating units from each animal in control conditions (black circles, n = 7) or 918 various methods of dopamine depletionbilateral 6OHDA (green diamond, n = 9), unilateral 919 6OHDA (blue triangle, n = 5), or systemic reserpine (orange square, n = 7). Lines indicate 920 mean. a. Delta (0.5-4 Hz) oscillations detected using both PSD peak and low phase shift 921 criteria. ANOVA: p = 5.206*10 -5 ; bilateral: p = 9.506*10 -5 ; unilateral: p = 0.00172; reserpine: p = 922 5.908*10 -5 , Dunnett's post-hoc test. b. Same as a, but using only the spectral power criterion. 923 ANOVA: p = 4.668*10 -4 ; bilateral: p = 5.645*10 -4 ; unilateral: p = 0.00601; reserpine: p = 924 5.794*10 -4 . c-d. Same as a-b, but for beta (7 -35 Hz) oscillations. With phase shift, ANOVA: p 925 ECoG. Lines between the bottom two panels illustrate M1 exhibiting peaks at a consistent time 972 for each animal in control (black circle, n = 7) and bilaterally dopamine depleted with M1 lesion 990 (dark green diamond, n = 3) conditions (p = 1.1478*10 -6 , two-sample t-test.) Fraction of oscillating units in SNr for each control animal (black circle, n = 7) and in the intact 996 hemisphere of unilaterally depleted animals (dark blue triangle, n = 4). The difference between 997 these conditions is not significant at the α = 0.05 level (p = 0.1138, two-sample t-test).