Distinctive Effects of D1 and D2 Receptor Agonists on Cortico-Basal Ganglia Oscillations in a Rodent Model of L-DOPA-Induced Dyskinesia

L-DOPA-induced dyskinesia (LID) in Parkinson’s disease has been linked to oscillatory neuronal activities in the cortico-basal ganglia network. We set out to examine the pattern of cortico-basal ganglia oscillations induced by selective agonists of D1 and D2 receptors in a rat model of LID. Local field potentials were recorded in freely moving rats using large-scale electrodes targeting three motor cortical regions, dorsomedial and dorsolateral striatum, external globus pallidus, and substantial nigra pars reticulata. Abnormal involuntary movements were elicited by the D1 agonist SKF82958 or the D2 agonist sumanirole, while overall motor activity was quantified using video analysis (DeepLabCut). Both SKF82958 and sumanirole induced dyskinesia, although with significant differences in temporal course, overall severity, and body distribution. The D1 agonist induced prominent narrowband oscillations in the high gamma range (70–110 Hz) in all recorded structures except for the nigra reticulata. Additionally, the D1 agonist induced strong functional connectivity between the recorded structures and the phase analysis revealed that the primary motor cortex (forelimb area) was leading a supplementary motor area and striatum. Following treatment with the D2 agonist, narrowband gamma oscillations were detected only in forelimb motor cortex and dorsolateral striatum, while prominent oscillations in the theta band occurred in the globus pallidus and nigra reticulata. Our results reveal that the dyskinetic effects of D1 and D2 receptor agonists are associated with distinct patterns of cortico-basal ganglia oscillations, suggesting a recruitment of partially distinct networks. Supplementary Information The online version contains supplementary material available at 10.1007/s13311-022-01309-5.


Introduction
Parkinson's disease (PD) is a progressive neurodegenerative disorder characterised by a severe loss of dopaminergic innervation to cortico-basal ganglia motor networks. This leads to an emergence of oscillatory neuronal activities that are believed to disrupt the dynamic processing of movement-related information and thus generate motor symptoms [1]. Local field potential (LFP) recordings during parkinsonian hypokinetic states demonstrate a synchronised oscillatory activity in the beta range  in both PD patients [2][3][4] and animal models of PD [5][6][7][8]. The exaggerated beta-band activity is suppressed by dopaminergic drugs as they improve parkinsonian motor symptoms [1] (for an extensive review of oscillatory activities in PD, see [9,10]). The dopamine precursor L-DOPA is the most efficacious symptomatic medication for PD, although it leads to a development of motor fluctuations and dyskinesia in the majority of patients (reviewed in [11]). These complications are associated with the emergence of new types of oscillatory neuronal activities in the basal ganglia. During the expression of L-DOPA-induced dyskinesia (LID), low-frequency oscillations in the theta range (5-10 Hz) have been detected in deep basal ganglia nuclei in both PD patients [12][13][14] and animal models of LID [15][16][17]. Moreover, recordings from the motor cortex have revealed a strong association between the expression of involuntary movements and characteristic narrowband gamma (NBG) oscillations, which have been detected in both the rat model of LID [18][19][20] and in PD patients [21]. Because of their strong association with dyskinesia, NBG oscillations have been proposed to serve as a feedback signal for adaptive deep brain stimulation (DBS) [22]. In general terms, NBG oscillations of the LFP are thought to reflect rhythmic synchronisations of transmembrane currents among populations of neurons in the recorded brain region (reviewed in [23]). However, little is known about the mechanisms through which L-DOPA elicits NBG oscillations specifically in dyskinetic subjects.
In both human PD and dopamine-denervated animals, L-DOPA exerts its motor effects by stimulating D1 and D2 receptors, which are highly expressed in the striatum but also present in other critical nodes of the cortico-basal ganglia network. In the rat motor cortex, the distribution of D1 and D2 receptors exhibits a rostral-caudal gradient of decreasing density, correlating with the distribution of dopaminergic terminals [24,25]. In the striatum, the expression of D1 and D2 receptors is segregated between "direct pathway" and "indirect pathway" spiny projection neurons (dSPNs and iSPNs, respectively), which influence the basal ganglia output via different routes (reviewed in [26]). Thus far, the role of these main dopamine receptor systems in the expression of LID-related oscillations has only been investigated in the primary motor cortex (forelimb area) and only upon acute administration of D1 or D2 receptor agonists, which were both reported capable of inducing cortical high gamma activity [19]. We here present a comprehensive investigation of oscillatory responses to D1 vs. D2 receptor agonists utilising a technology that enables large-scale multi-structure recordings in awake behaving animals [27]. Applying this technology to a validated rat model of LID [28], we analyse LFP oscillations in both NBG and theta frequency bands from seven critical nodes of the cortico-basal ganglia motor network. Our results reveal prominent differences in the oscillatory patterns mediated by D1 vs. D2 receptors in dyskinetic animals, suggesting a recruitment of partially distinct networks.

Animals
Experiments were performed in adult female Sprague-Dawley rats, as in our previous studies on the subject [15,16,18]. Animals (purchased from Janvier Labs, France) weighed 250 g on arrival and were housed in standard cages under controlled temperature (22 °C) and humidity (50%) laboratory conditions on a 12:12 h light/dark cycle (lights on at 6:00 a.m.). Food and water were available ad libitum. All procedures were approved by the Malmö-Lund Ethical Committee on Animal Research.

Study Design
Rats sustained unilateral 6-hydroxydopamine (6-OHDA) lesions and 3 weeks after the lesion, they were primed with L-DOPA (6 mg/kg s.c.) daily for 1 week and then every second day for 2 weeks before electrode implantation surgery (Fig. 1A). Starting 1 week from electrode implantation, all animals received treatment with L-DOPA (n = 9), the D2 receptor agonist sumanirole (2 mg/kg s.c., n = 8), and the D1 receptor agonist SKF82958 (0.05 mg/kg s.c., n = 9), according to the order shown in Fig. 1A, with 2-3 injections per week for a total of 6-8 weeks. Animals were recorded on vehicle treatment with saline at different time points throughout the drug treatment period (veh; vehicle; n = 9). Each treatment was administered every other day for 1-2 weeks and electrophysiological recordings and dyskinesia rating sessions were carried out on the days of drug challenges (Fig. 1A). However, when introducing a new treatment, a drug-free break of minimum three days was introduced between treatments and the new drug was administered at least twice before recordings were recommenced to ensure a stable level of dyskinesia. All animal were perfused and CT scans and immunohistochemical analysis were performed to identify the location of the tip of each wire ( Fig. 1A-E).

Drug Treatments
Dose-response curves were performed in a separate group of non-implanted rats prior to commencing the physiological recording study (see Supplemental Methods). The doses chosen were well tolerated (Supplemental Figs. S1 and S2) and in the range of those previously reported to produce robust behavioural effects in rodents through the indicated targets [29][30][31][32][33][34]. L-DOPA (levodopa methyl ester hydrochloride, 6 mg/kg; Sigma Aldrich AB, Sweden), sumanirole (SUM; sumanirole maleate, 2.0 mg/kg; Tocris, UK), and SKF82958 (SKF-82958 hydrobromide, 0.05 mg/kg; Tocris, UK) were given subcutaneously (s.c.). L-DOPA was co-administered with the peripheral DOPA decarboxylase inhibitor benserazide hydrochloride (12 mg/kg s.c., Sigma Aldrich AB, Sweden). All drugs used in this study were dissolved in physiological saline solution and administered in a volume of 1 ml/kg.
Three weeks after the lesion, animals with severe (> 85%) unilateral dopamine denervation were selected using the cylinder test of forelimb use asymmetry [35,[37][38][39]. Briefly, rats were placed individually in a glass cylinder (25 cm in diameter and 40 cm in height) and video filmed for 3-5 min. The total number of supporting wall contacts performed independently with the left and right forepaw was counted Experimental procedure and validation of wire placement. A Schematic representation of the study design. Rats sustained 6-OHDA lesions in the right MFB and were primed with L-DOPA (6 mg/kg s.c.) before implantation surgery. When recovered, all animals received treatment with L-DOPA, the D2R agonist sumanirole (SUM; 2 mg/kg s.c.), the D1R agonist SKF82958 (0.05 mg/kg s.c.), and vehicle (veh) with 2-3 injections per week for a total of 6-8 weeks. A drug-free break of three days was introduced between treatments. When introducing a new dyskinetic treatment, the drug was administered at least twice before recordings were recommenced to ensure a stable level of dyskinesia. Following perfusion, CT scans and immunohistochemical analysis were performed to identify the location of the tip of each wire. B Representative tyrosine hydroxy-lase (TH) immunostaining of dopaminergic neurons in a horizontal section of one animal. C Dorsal view of the rat skull with a constellation of wire targets. (Adapted from Paxinos and Watson [45]). D Horizontal (top) and sagittal (bottom) view of the location of the tip of each wire identified from CT scans. E Comparison between immunohistochemical analysis of CV stained brain sections and CT segmentation overlaid the corresponding atlas from Paxinos and Watson [45] for SNr of one example brain. Arrows indicate one electrode tip identified by each method. Abbreviations: computed tomography (CT), cresyl violet (CV), rostral forelimb area (RFA), primary motor cortex forelimb area (M1FL), primary motor cortex trunk area (M1Tr), dorsomedial striatum (DMS), dorsolateral striatum (DLS), globus pallidus pars externa (GPe), and substantia nigra pars reticulata (SNr) offline. Then, counts from the forepaw contralateral to the lesion were expressed as a percentage of the total number of wall contacts. Starting 3 weeks after the lesion, rats with less than 25% contralateral paw usage were primed with L-DOPA (Fig. 1A). Recording electrodes were implanted in animals that developed moderate to high levels of dyskinetic behaviours (cut-off criterion: AIM severity grade ≥ 2 on at least two AIM subtypes on most monitoring periods, as in [40]). In total, 11 animals were selected for electrode implantations surgery using these criteria.

Abnormal Involuntary Movements (AIMs) Ratings
Dyskinesia was assessed online using both the basic AIMs scale [36] and an amplitude scale [35]. The basic scale evaluates AIM severity based on the percentage of observation time during which a dyskinetic feature is present, considering three subtypes of dyskinetic movements (axial, forelimb, and orolingual AIMs) [36]. The basic scale was combined with the AIMs amplitude scale, which evaluates the deviation of a body part from its natural resting position and the number of muscle groups visibly engaged in the dyskinetic movement [35]. Thus, all AIM subtypes were scored simultaneously with respect to their severity and amplitude following drug injection for monitoring periods of 1 min every 5 min for the first 20 min and then every 10 min for the rest of the testing session. Each of the AIM subtypes received a severity score and an amplitude score in a scale from 0 to 4 on each monitoring period. A composite AIM score was then produced by multiplying the severity score and amplitude score for each AIM subtype on each monitoring period and all AIM scores were summed to obtain a global measure of dyskinesia for each testing session (global AIM score).

Recording Electrodes
A recording implant of 94 insulated tungsten wires (33 µm diameter; California Fine Wire Company, CA, USA) was built as previously described [27]. The microwires were arranged into 14 electrode bundles (7 per hemisphere) with 250 µm spacing between wires in each horizontal dimension and a wire length corresponding to the depth of each recording target (see Table 1 for coordinates and main references) (Fig. 1C). A 125-µm silver wire (Advent Research Materials Ltd, UK) was used for ground connection. All wires from each hemisphere were connected to a customdesigned printed circuit board with silver conductive paint (Electrolube, UK) and secured with UV-curable adhesive (Dymax, CT, USA).

Implantation Surgery
Electrode implantation surgeries were performed at least 6 weeks after 6-OHDA lesions (Fig. 1A). Animals were anaesthetised with fentanyl/medetomidine (0.21/0.21 mg/ kg, i.p.; Apoteket AB, Sweden) and fixed in a stereotaxic frame (David Kopf Instruments, CA, USA) in a flat skull position. Holes were drilled in the frontal and parietal skull bone over the target recording sites followed by durotomy to accommodate a bilateral implantation of the microwire arrays. Silver wire was attached to five screws in the frontal and occipital skull bone for ground connection. The implant was anchored to the skull screws with dental acrylic cement (Kerr, CA, USA) also covering the ground wires for electrical insulation. After surgery, the anaesthesia was reversed by atipamezole hydrochloride (0.5 mg/kg, i.p.; Apoteket AB, Sweden), and buprenorphine (0.05 mg/kg, s.c.; Apoteket AB, Sweden) was administered as postoperative analgesic. The animals were allowed to recover for a minimum of 1 week before recording sessions commenced. During this period, they received extra postoperative care until fully recovered.

Electrophysiological Recording Sessions
For one animal, a flat signal was obtained in all recordings and it was excluded from the electrophysiological analysis. In total, 79 recordings from 10 unilaterally 6-OHDAlesioned dyskinetic rats were included in the study, of  [71] which 7 animals completed the entire series of treatments. During recording sessions, animals were placed individually in a transparent plastic cylinder (diameter 55 cm, height 40 cm) and their behaviour was recorded with digital video from the side and above (25 frames per second, Genie HM640 GigE camera; Teledyne DALSA, Canada) in parallel with electrophysiological signals (see the "Signal Acquisition" section). Synchronisation between video frames and electrophysiological signals was performed using a Master-8 pulse generator (A.M.P.I, Israel). All recording sessions followed the same timeline. After being placed in the cylinder, animals were allowed to habituate for 10 min before they were briefly anaesthetised with isoflurane (5%) to connect custom-build adapters, Intan recording headstages (Intan Technologies, LLC, California, CA, USA), and recording cables to the implant. After recovering from the anaesthesia, animals were recorded for 30 min to establish baseline conditions and then injected with vehicle. After 20 min of recording, animals received an injection of either L-DOPA, SKF82958, sumanirole, or vehicle followed by a further 3 h of recording before ending the experiment. For details on the pharmacological paradigm see the "Drug Treatments" section.

Signal Acquisition
Wideband neuronal activity was recorded with the Open Ephys acquisition system [41] using four Intan RHD2132 amplifier chips with on-board AD converters (Intan Technologies, LLC, California, CA, USA). The wideband signals were bandpass filtered between 1 and 7600 Hz and digitised and recorded at 30 kHz. Then, LFPs were extracted off line, low-pass filtered (8th order Butterworth at 500 Hz) and downsampled to 2000 Hz.

Time-Frequency Analysis of LFP Power
To emphasise the local sources of the measured electrical potential, bipolar LFP time series from all unique pairs of wires within the same structure were computed offline for each LFP recording using custom code written in MATLAB (MathWorks). For each of these time series, a spectrogram (i.e. time series of power spectral densities (PSDs)) was estimated over the 0-250 Hz frequency range with Welch's method (8 s Hann window with 50% overlap). To separate oscillatory components in the power spectrum from the arrhythmic, so-called fractal (2) These mixed times series (composed of both fractal and oscillatory signals) were used as input for irregular-resampling auto-spectral analysis (IRASA). By irregularly resampling the time series multiple times, the fractal components could be separated from the mixed time series (left). The PSD was the normalised to the fractal component to construct a power spectrum measure that emphasises truly rhythmic activity (right). (3) To detect oscillations in each frequency band of interest, we fitted a parametric model to the spectral peak. Examples of LFP power spectra (grey) from 8 s of data are shown and green or red lines show the corresponding fitted function y(f). Functions with parameter values and goodness-of-fit within given limits were counted as successful peak detections and marked in green. (4) Example of a spectrogram where each green point marks an 8-s bin with successful peak detection. Hatched pink square marks the segment of spectrum fitted to only detect the frequency band of interest. Examples for each treatment can be found in Supplemental Fig. S3 components, we used the irregular-resampling auto-spectral analysis (Wen and Liu [42]) (Fig. 2). By irregularly resampling the time series multiple times and then normalising the spectrum to the fractal component, it was possible to construct a power spectrum measure that emphasises truly rhythmic activity: where S(f) and S fractal (f) have the dimension power per frequency (i.e.V 2 /Hz), and S dB(fractal) is expressed in the dimensionless unit dB fractal . As a final step, spectrograms from each individual local bipolar LFP time series were averaged to create one representative spectrogram from each structure. Channels with exceptional noise levels were excluded based on visual inspection. Also, power line noise (50 ± 2 Hz and harmonics at 100 ± 1 Hz, 150 ± 1 Hz and 200 ± 1) was removed from the PSDs.

Detection of Oscillatory Activity
To detect oscillations in the theta, beta, and NBG bands, we defined a parametric model and fitted it to the spectral peak (Fig. 2). The parameters A (peak height), B (peak frequency), C (peak width), D (inclination of flat background), and E (offset of flat background) were estimated such that the model y(f) fit the spectrum optimally in the least-squares sense. The reliability of the estimated model was assessed by the goodness of fit (R2). Since the model only included one peak, different segments of the spectrum were fitted depending on the frequency band of interest: [1 15] Hz for theta, [8 48] Hz for beta, and [40 120] Hz for NBG (Fig. 2). Limit values for the fitted parameters and R2 were defined through visual inspection, and if the parameters for a given spectrum fell inside those limits, a positive detection was declared ( Fig. 2; see Supplemental Table S1 for typical limit values). This procedure resulted in a possible detection each 4 s. The detection rate was then calculated as the average of the binary outcome for each treatment within the structures of interest. If sufficient detections were made (> 5%), further analysis of the oscillation intensity (the height of the peak) and oscillation frequency was performed.

Phase Synchrony and Functional Connectivity
Phase coherence was investigated on simultaneously recorded NBG signals from all structures in which these oscillations were detected during peak dyskinesia, and for all treatments. Signals were initially bandpass-filtered ± 5 Hz around the median frequency. Then, the Hilbert transform of each signal was computed with the MATLAB function hilbert (MathWorks). The instantaneous phase angle (t) was then obtained by calculating (t) = argz(t) using the MATLAB function angle (MathWorks). The mean phase difference between a given pair of wires i and j was calculated as Δ ij = ⟨ i (t) − j (t)⟩ , where ⟨⟩ denotes the circular mean (circ_mean, MATLAB CircStat toolbox; [43]). Similarly, the resultant vector length r ij was obtained using the function circ_r of the same toolbox.
To obtain a measure of phase coherence between two structures, r ij for all unique pairs with wire i in the first structure and wire j in the second structure were averaged. Similarly, the mean phase difference between two structures was obtained by averaging all unique Δ ij . However, to obtain a good quality phase estimate, we only included Δ ij for pairs with a mean resultant vector length > 0.5, assuming that a higher magnitude reflects a higher signal to noise ratio of the source of interest and thus a better-quality phase estimate.

Video Tracking Analysis
The position of each rat was tracked from the top camera video recordings with the open-source tool DeepLabCut [44]. Briefly, a total of 1221 frames picked randomly from 23 videos with different brightness conditions were used to label different body parts of the rat in order to train a deep learning network over 1,030,000 iterations. The network was then fed with all top camera videos from this study to track the chosen body parts.
Next, x and y coordinates of 7 body parts, that is, nose tip, left and right ear, head top (middle point on the implant between both medial headstages), body center (point on the midline running on the back between the thorax and abdomen), base of tail, and tip of tail, were extracted using MATLAB (MathWorks), considering only a detection likelihood of > 0.9. The centroid of the animal was calculated by averaging the x and y coordinates of all 7 body parts and the front of the animal by averaging x and y coordinates from nose, left and right ear, and head top. The direction of the animal was defined as the direction of the line between the front of the animal and the base of the tail.
To reveal the egocentric rotational behaviour of the rat, the number of full 360° turns was counted as follows: The direction measurements were binned into 8 segments. If the direction measurements passed through all 8 sectors in the correct order, it was counted as a full turn. Contralateral and ipsilateral turns were counted separately. An incomplete rotation reset the counting of sector transitions and the last angle in the discarded rotation was used as the first angle of the next count. The start and end frames of each full turn were used to determine the distance travelled during the turn and the duration of the turn.
Overall motor activity was computed by dividing the recording chamber into 1 × 1 cm squares and counting the number of squares visited in 1 min recording bins with reference to the body centroid. The number of squares visited was then divided by the total number of squares, and the corresponding value was represented as "fraction of total area visited". In addition, the speed of overall motor activity was calculated based on the translation of the body centroid during a 1 s window. Overall motor activity, overall motions, and open field motions are used interchangeably to describe motor activity in the open field arena.

Computed Tomography (CT) Scanning of the Implanted Microwire Arrays
The locations of recording wires were identified using computed tomography (CT) scanning of the heads with the microwire arrays still implanted. CT scans were performed on a Mediso NanoScan PET/CT scanner (Mediso, Hungary). To minimise CT metal-artefacts, the animal head was positioned such that the wires were perpendicular to the photon beam. This was achieved by placing the head on its side inside the tunnel such that the tunnel axis was aligned with the DV direction. The CT was performed as a helical scan with 0.30 pitch, 700 projections with 300 ms exposure time, 70 kV tube voltage, and 72 µA tube current. Images were reconstructed with thin-slice Ram-Lak filter, to give a 30 µm side, and 34 µm slice thickness (voxel size). The tip of each wire was manually identified in the resulting images. Bregma and lambda were used as landmarks to register the images to an anatomical atlas of the rat brain [45], followed by a manual calibration to visually optimise the alignment between atlas and scan. The atlas coordinates of the wire tips were calculated from the voxel coordinates using the resulting affine transformation. Finally, each wire was assigned an appropriate anatomical label based on the location of the wire tip in the atlas ( Fig. 1D-E).

Histological Verifications
Animals were anaesthetised with a lethal dose of sodium pentobarbital (100 mg/kg i.p., Apoteksbolaget AB, Sweden) and transcardially perfused with 0.9% saline followed by 4% paraformaldehyde before decapitation. Heads were kept in paraformaldehyde at 4 °C for post-fixation. After CT scanning of the head for electrode placement verification (see the "Computed Tomography (CT) Scanning of the Implanted Microwire Arrays" section), the brains were extracted and transferred to 25% sucrose solution in 0.1 M phosphate buffer (PB) at 4 °C until sinking for cryoprotection (at least 1 day). Coronal sections of 30 µm were cut serially through the entire brain using a microtome and stored at − 20 °C in non-freezing buffer solution until further processing for immunohistochemistry. The extent of dopamine denervation was verified in each animal by immunohistochemical staining for tyrosine hydroxylase (TH) (see Supplemental Methods). The placement of electrodes was evaluated in cresyl violet-stained sections by microscopic analysis and compared to CT scans of the implants to verify the locations of recording wires (Fig. 1E). Wires that did not target the structures of interest were excluded from further analysis.

Statistical Analysis
Statistical analysis of behavioural data was performed using Prism 9 (GraphPad Software). Comparisons of treatment effects between groups, on single sessions or time points, were done with a mixed-effects model with Geisser-Greenhouse correction, followed by post hoc Tukey's multiple comparisons test. All behavioural data compared between groups on single sessions or time points are presented as box plot and median, with whiskers annotating minimum and maximum values.
Detection rates of LFP frequency bands were compared with a generalised linear mixed-effects model (R glmer function in lme4 package) with binomial distribution and 4 factors (treatment, structure, recording, and animal). The factors recording and animal were considered random with recording nested in animal. Hemisphere was included as factor when appropriate. Bias adjustment was applied to account for the independent sources of random variation estimated by the model. Scheffé's multiple comparisons post hoc test was chosen and results are reported as odds ratios (OR) (details on multiple comparisons can be found in Supplemental statistical materials-Detections).
Comparisons of peak height (absolute power) and peak frequency were done with linear mixed-effects models (R lmer function in lme4 package) with log-normal distribution when appropriate and factors treatment, structure, recording, and animal. The factors recording and animal were considered random with recording nested in animal. Scheffé's multiple comparisons post hoc test was chosen and results are reported as estimates of differences or ratios (details on multiple comparisons can be found in Supplemental statistical materials -Peak height or Supplemental statistical materials -Peak frequency).
Comparisons of instantaneous phase differences were done with a MANOVA approach based on trigonometric functions of the circular data [46] using RStudio with factors treatment and structure. In addition, a one sample test for the mean angle was applied (circ_mtest in MATLAB Circ-Stat toolbox [43]). Scheffé's multiple comparisons post hoc test was chosen and results are reported as estimates of differences. Comparisons of phase coherence were done with Wilcoxon signed-rank test (MATLAB signrank function) (details on multiple comparisons can be found in Supplemental statistical materials -Phase or Supplemental statistical materials-Functional connectivity).
All electrophysiological parameters are given in mean ± SEM unless otherwise stated. Relations between variables were examined using the Spearman's Rho. The level of statistical significance was set at α = 0.05 and p-values for detection rates, peak height, and peak frequency were adjusted to obtain multiplicity-adjusted results. To account for multiple testing in phase analysis, the Bonferroni correction was applied setting the significance at α/n.

Patterns of Dyskinesia and Rotational Behaviour Induced by D1 vs. D2 Receptor Stimulation
All 6-OHDA-lesioned rats in this study were primed with daily L-DOPA treatment at the dose of 6 mg/kg/day, which corresponds to a therapeutic-like dose for this PD model [47][48][49]. After implanting the recording electrodes, all animals were treated sequentially with L-DOPA, followed by 2.0 mg/kg sumanirole (D2 receptor agonist), and then 0.05 mg/kg SKF82958 (D1 receptor agonist) over a period * * * * * * of several weeks (Fig. 1A). The doses of sumanirole and SKF82958 were carefully selected based on pilot experiments in non-implanted rats (Supplemental Fig. S2) as well as previous literature [29][30][31][32][33][34]. These doses represent a middose range, where off target effects of the selected compounds can be assumed to be negligible. Dyskinetic behaviours were induced by all the compounds tested, although with differences in time course and overall severity ( Fig. 3A-C). L-DOPA induced dyskinesia in all animals, reaching peak severity at 50-60 min post injection (Fig. 3A), with AIM scores comparable to those previously reported from non-implanted rats [48,49]. SKF82958 resulted in an immediate induction of dyskinesia yielding a severity peak at 20-30 min post injection, though with AIM scores similar to those recorded during the severity peak of L-DOPA (Fig. 3A, C). However, the time course was significantly shorter compared with both L-DOPA and sumanirole, and AIM scores showed an abrupt decrease down to baseline between 70 and 100 min post injection (Fig. 3A). After treatment with sumanirole, both total and peak AIM scores were significantly lower compared to both L-DOPA and SKF82958 (Fig. 3B-C; Total: − 41% vs L-DOPA, *p = 0.008; − 30%, vs. SKF82958, #p = 0.004; Peak: − 54% vs. L-DOPA, *p = 0.001; − 54% vs SKF82958 #p < 0.001). Similar to L-DOPA, peak AIM severity with sumanirole was seen at 50-60 min post injection (Fig. 3A). In order to compare the treatments during functionally analogous time windows, we focused our analyses on the periods corresponding to peak AIM severity, i.e., 20-60 min post injection for SKF82958 and 40-80 min for both L-DOPA and sumanirole (see coloured bars in Fig. 3A).
Analysis of the individual AIM subtypes revealed that axial, limb, and orolingual scores were all significantly lower after sumanirole treatment compared with L-DOPA and SKF82958 (Supplemental Fig. S4A-C, E-G). Significant differences between the agonists and L-DOPA were detected upon expressing each AIM subtype scores as a percentage of the total scores recorded from each animal during the test session. Compared with L-DOPA, sumanirole treatment induced a higher relative expression of axial AIMs and lower expression of orolingual AIMs (Supplemental Fig. S4D, H).
Investigating rotational behaviour, we found a difference in the time course between treatments (Fig. 3D). Indeed, both dopamine agonists induced more contralateral turns compared to L-DOPA ( Taken together, these results show that dyskinetic behaviours developed upon pharmacological stimulation of both D1 and D2 receptors in L-DOPA-primed animals, although dyskinesia severity was milder upon treatment with the D2 agonist compared to both L-DOPA and SKF82958. In addition, both dopamine agonists induced more contralateral turns compared to L-DOPA, these turns being faster upon D1 receptor stimulation.

Dyskinetic Behaviours Concur with Network-Wide Changes in LFP Activity Patterns
In parallel with the behavioural assessment, LFP recordings were performed simultaneously in 7 nodes of the corticobasal ganglia network of both the intact and the lesioned hemisphere. Indeed, we recorded from the rostral forelimb area (RFA), two primary motor cortical areas suggested to be involved in forelimb (M1FL) and trunk (M1Tr) movements, the dorsolateral (DLS) and dorsomedial (DMS) striatum, and the primary targets of striatal projection systems, that is the globus pallidus pars externa (GPe) and the substantia nigra pars reticulata (SNr) ( Table 1). An example multi-site recording for each treatment is illustrated in Fig. 4B with the time course of the corresponding AIM scores and speed of open field motions presented in the bottom panel. We noted a close association between the expression of dyskinetic behaviours and the occurrence of NBG LFP oscillations  in several of the recorded structures in the lesioned hemisphere. Spectral contents were in fact altered in several frequency bands during peak-of-dose dyskinesia (Fig. 4B). In particular, beta oscillations (12-35 Hz) were generally decreased compared to vehicle and baseline conditions (see Fig. 4A and the interval preceding drug injection in Fig. 4B). The occurrence of beta oscillations in PD has been covered extensively elsewhere [3-7, 12, 19, 50-53], and has furthermore shown little association to dyskinesia [19]. Therefore, beta oscillations will not be investigated further in this study. In addition to NBG oscillations, several of the recorded structures exhibited a prominent increase in theta oscillations (5-10 Hz). A summary of the LFP power spectra during peak-of-dose dyskinesia is shown in Fig. 4C, which indicates conspicuous increases in both the NBG (70-110 Hz) and the theta range (5-10 Hz) in most of the recorded structures.
In the following, we will therefore focus our analyses on these frequency bands and compare the effects of different treatments during peak AIM severity, i.e. 20-60 min post injection for SKF82958, 40-80 min for L-DOPA, sumanirole, and vehicle (cf. coloured bars in Fig. 3D).

Differential Effects of D1 and D2 Receptor Agonists on Narrowband Gamma (NBG) Oscillations (70-110 Hz)
Following treatment with L-DOPA and the dopamine agonists, NBG oscillations appeared in several of the recorded structures of the lesioned hemisphere coinciding temporally with the AIMs (Fig. 4B-C and Supplemental Fig. S5A), which is in line with previous findings in animal models of LID and reports from dyskinetic PD patients [12,16,18,19,21]. NBG oscillations were not detected in the intact hemisphere, nor were they detected after vehicle treatment in either hemisphere ( Fig. 5A and Supplemental Fig. S6). The below report is therefore focused on a comparison between the three pharmacological treatments in the lesioned hemisphere.
To further emphasise the relation between LFP oscillations in the NBG range and dyskinesia, we calculated the correlation of NBG power against manually scored AIMs. Indeed, NBG oscillations were strongly correlated with AIM scores in the majority of structures recorded, being strongest for M1FL (Rho = 0.680, p < 0.001, Supplemental  Fig. S8A). Significant positive correlations were moreover found between NBG power and each individual AIM subtypes (see Supplemental Table S2).
Because of the correlation found between NBG oscillations and dyskinesia severity, we examined this oscillatory activity at time points when the dopamine receptor agonists were producing equivalent AIM scores (i.e. 70-80 min after the administration of SKF82958/sumanirole). At these time points too, NBG detection rates in subcortical structures were significantly lower after sumanirole compared to SKF82958 administration (Supplemental Fig. S7A, B) and relative differences as to NBG power were maintained in most structures (Supplemental Fig. S7C). This suggests that the stronger action of D1 vs. D2 receptor stimulation in driving NBG oscillations is not contingent on an induction of more severe dyskinesia.

D2 Receptor Stimulation Elicits Theta Oscillations (5-10 Hz) in the Deep Basal Ganglia Nuclei
Theta oscillations were present in all recorded structures in both the intact and lesioned hemispheres and had a mean frequency of 6.8 ± 1.2 Hz. The detection rate of theta oscillations in the lesioned hemisphere across structures was increased after treatment with vehicle, L-DOPA, and SKF82958 compared to the intact hemisphere, whereas sumanirole induced a bilateral increase in the detections of theta oscillations (Mixed-effects model: χ 2 (treatment) 3 = 56.7, p < 0.001; χ 2 (hemisphere) 1 = 733, p < 0.001; χ 2 (interaction) 3 = 364, p < 0.001; OR veh intact/veh lesioned = 0.63 ± 0.02, p < 0.001; OR L-DOPA intact/L-DOPA lesioned = 0.80 ± 0.01, p < 0.001; OR SKF82958 intact/SKF82958 lesioned = 0.84 ± 0.01, p < 0.001; OR SUM intact/SUM lesioned = 1.1 ± 0.02, p = 1.0) (Fig. 6A-B). Furthermore, the detection rate of theta oscillations was constant during the time course of the dyskinetic behaviour with only small changes in the early and late phase of dyskinesia ( Fig. 4B and Supplemental Fig. S5B). When investigating the detection of theta oscillations within the cortical structures, the detection rate was increased in the M1FL area of the lesioned hemisphere by all treatments compared to vehicle ( Fig. 6A; M1FL: OR veh/L-DOPA = 0.21 ± 0.03, &p < 0.001; OR veh/SKF82958 = 0.36 ± 0.05, &p < 0.001; OR veh/SUM = 0.26 ± 0.04, &p < 0.001). In addition, treatment with L-DOPA increased the detection rate of theta oscillations bilaterally in both RFA and M1FL as well as in the striatum compared to all other treatments (Fig. 6A-B and Supplemental Table S4;  Thus, there was an increased detection rate of theta oscillations after treatment with sumanirole in the deeper nuclei (GPe and SNr) compared to the other treatments, but the detection rate in the GPe was also increased after treatment with SKF82958 and L-DOPA compared to vehicle.
To explore the relation between theta oscillations and dyskinesia, we calculated the correlation of theta power against global AIM scores. Weak, though significant, correlations were found in RFA, M1FL, and striatum (DMS and DLS) (Supplemental Fig. S8B). Noteworthy, a robust correlation between dyskinesia severity and theta power was found in the GPe of the lesioned hemisphere, where both the global AIM scores (Supplemental Fig. S8B; Rho = 0.338, p < 0.001) and the individual AIM subtypes (Supplemental Table S3) were significantly correlated with the power of theta oscillations.
Since increases in theta activity had also occurred on the intact side, we moreover analysed the correlations between theta power and an index of overall motions in the recording chamber (the fraction of total area visited). We focused this analysis on the GPe and SNr, where the largest increases in theta detection rate and power had been detected after treatment (see Fig. 6A-D). In both GPe and SNr, there was a positive correlation between the fraction of area visited and theta power in both the intact and lesioned hemisphere across treatments (Supplemental Fig. S9A, B). However, the relationship between overall motions and theta power became negative in the lesioned GPe specifically upon treatment with sumanirole (Supplemental Fig. S9A; Rho = − 0.265, p = 0.002). This indicates that theta power measured in the dopamine-denervated GPe after sumanirole treatment was specifically correlated with the expression of dyskinesia, not with an increase in motor activity. In contrast, treatment-induced bilateral theta-increases in the SNr were not directly related to the presence of dyskinesia (cf. Supplemental Fig. S8B) but seemed to reflect an overall increase in motions (cf. Supplemental Fig. S9B).

NBG Phase Synchrony Reveals a Stronger Effect of D1 Stimulation on Cortico-Subcortical Coupling
To identify potential drivers of NBG oscillations in the corticobasal ganglia network, we examined phase relationships between different structures during peak dyskinesia, as induced by L-DOPA, sumanirole, and SKF82958 (Fig. 7A). This revealed that the phase differences between NBG oscillations in different structures were not random, as would be expected if there was no functional coupling between the structures. Rather, we found a strong bias towards certain phase difference values. Using the strength of this bias as a measure of functional connectivity, we found a stronger functional coupling between structures in the lesioned hemisphere compared to the intact hemisphere (Supplemental Fig. S10). Notably, a tighter functional coupling was seen within cortex and between cortex and striatum in the lesioned hemisphere upon treatment with SKF82958 compared to sumanirole  Table S6). While a tight functional coupling was found between striatum and GPe following all treatments, the GPe also showed strong functional connectivity with M1Tr and RFA after treatment with SKF82958/L-DOPA and L-DOPA, respectively (Fig. 7B-D and Supplemental Table S6).
Generally, NBG phase differences between structures were concentrated around 0°. However, several structure pairs showed a non-zero phase difference. As shown in Fig. 7A, M1FL was leading RFA and striatum for all treatments ( Fig. 7A and Supplemental Table S7), with a differential effect depending on treatment. NBG oscillations in RFA were moreover lagging GPe for both SKF82958 and  Table S7). For the remaining structures, NBG oscillations were phase synchronised. Altogether, our data show that all recorded structures are strongly coupled with small phase differences < 30°, corresponding to delays < 1 ms. Also, the fact that M1FL leads RFA and striatum suggests that it might be an important driver of NBG oscillations.

Discussion
Using in vivo multi-structure recordings in a rat model of LID, we have examined the pattern of cortico-basal ganglia oscillations induced by chronic treatment with potent and selective agonists of D1 and D2 receptors. The D1 agonist SKF82958 induced strong NBG oscillations in all recorded structures except for the SNr, whereas the D2 agonist sumanirole only induced NBG oscillations in forelimb motor cortical areas and dorsolateral striatum, and with significantly inferior power compared to SKF82958. In contrast, sumanirole had a stronger effect on the induction of theta oscillations, whose detection rate and power in the GPe and SNr were significantly larger than after SKF82958 treatment. To our knowledge, this is the first report on the effects of pharmacological D1 vs. D2 receptor stimulation on LFPs recorded simultaneously in many critical structures of the cortico-basal ganglia network. For the first time, we recorded simultaneously from motor cortical areas corresponding to forelimb (RFA and M1FL) or trunk (M1Tr) regions, from both associative and motor striatal regions (DMS and DLS, respectively), and from the primary targets of striatal projection systems that jointly control movement (GPe and SNr).

Effects of D1 and D2 Receptor Agonists on NBG Activity
Previous studies in rat models of LID have shown that the expression of dyskinesia induced by L-DOPA coincides with the appearance of NBG oscillations in the motor cortex and several other nodes of the cortico-basal ganglia-thalamic network [15,16,18,19,53,54]. Interestingly, NBG oscillations have been described in the motor cortex and subthalamic nucleus of dyskinetic PD patients [21], and have been proposed as a promising control signal for adaptive DBS [22]. Of note, this NBG oscillation was minimally affected by voluntary movements, while its presence was coinciding with dyskinesia [21], suggesting that it is distinct from movement-related broadband gamma synchronisation (reviewed in [55]). In the current study, sustained NBG oscillations  Hz) emerged in several cortico-basal ganglia structures in conjunction with dyskinetic behaviour induced by chronic L-DOPA treatment. Compared to other publications on the subject, rats were treated with a relatively low dose of L-DOPA (6 mg/kg s.c.) that is commonly used in behavioural studies [47][48][49]. In agreement with reports using a higher dose of L-DOPA (12-15 mg/kg) [15,16,[18][19][20] NBG oscillations were detected in the motor cortex of the lesioned hemisphere in all rats as well as in the striatum and GPe, whereas the detection rate in the SNr was limited. In contrast, a recent study did not find an increase in NBG oscillations in 6-OHDA-lesioned rats after treatment with 7 mg/kg L-DOPA in a 3-week priming paradigm, despite clear dyskinetic behaviour [56]. The discrepancy with our data could depend on differences in the methods of spectral analysis, where the use of IRASA in the current study (see the "Time-Frequency Analysis of LFP Power" section) allowed the isolation of true oscillatory phenomena [42].
The induction of NBG oscillations upon activation of D1 or D2 receptors has so far only been reported in specific motor cortical areas of the frontal lobe and upon acute administration of dopamine receptor agonists [18,19]. This is the first comprehensive investigation comparing the potency of chronically administered D1 and D2 receptor agonists at generating NBG oscillations in several cortical and subcortical areas. Moreover, whereas previous studies have used D2-class agonists with high activity on D3 receptors [19], we have here used the very selective D2 receptor agonist sumanirole. Stimulation of the D1 receptor elicited strong NBG oscillations not only in motor cortices but also in DMS, DLS, and GPe, that is, in all the recorded structures with the exception of SNr. Following treatment with the D2 agonist, NBG oscillations were elicited only in forelimb motor areas and DLS, but not in M1Tr, DMS, GPe, and SNr. In the latter structure, none of the tested treatments induced NBG oscillatory activity, indicating that the basal ganglia output via SNr is less prone to maintain these high-frequency oscillations. In the frequency domain, NBG oscillations were strongly modulated by the pharmacological treatments, where the mean frequency of the elicited NBG oscillations was highest after L-DOPA treatment and lowest after treatment with sumanirole. Similar modulations of broadband gamma oscillations have been reported in anaesthetised hemiparkinsonian rats, with higher peak frequencies after treatment with L-DOPA compared to the potent dopamine receptor agonist apomorphine [57].
In agreement with previous studies [18,19], we found strong NBG oscillations in forelimb motor cortical areas, which were more prominent than in all other structures recorded and correlated strongly with dyskinesia severity (AIM scores). Interestingly, applying a D1 receptor antagonist on the cortical surface has been found to stop both LID and cortical NBG oscillations in hemiparkinsonian animals [18], strongly indicating a direct cortical involvement. Thus, Halje et al. [18] proposed that loss of cortical dopaminergic innervations may be a key predisposing factor for dyskinesia and that the resulting sensitivity to dopamine makes the cortical circuits prone to network resonance in the NBG frequency range. As a novel contribution, our data indicate a particularly prominent role of the D1 receptor in inducing these cortical NBG oscillations, which is in line with the critical role of D1 receptor supersensitivity in LID [26,58]. Indeed, the D2 receptor agonist sumanirole was significantly less potent than SKF82958 at eliciting NBG oscillations. The superior potency of the D1 agonist might partly reflect the higher relative density of D1 receptors compared to D2 receptors in motor cortical regions [59], although further investigations are needed to clarify the specific cellular distribution of these two receptor classes in the cortex, and their changes following dopamine denervation. In the phase synchrony analysis, we found that NBG oscillations in the M1FL were leading motor cortical and striatal structures, with a stronger effect following treatment with the D1 receptor agonist (see further discussion below). Considering that M1FL showed the highest NBG detection rate, these findings point to this structure as a potential driver of these oscillations, at least for the circuits recorded in this study.
Most studies have so far focused on NBG oscillations in frontocortical motor areas, and it has remained unclear whether this rhythm could be induced in other cortical regions. For the first time, we here recorded NBG activity in a more caudal primary motor cortical region located in the parietal lobe (M1Tr), previously reported to be involved in trunk movements [60,61]. In this area, the induction of NBG oscillations was considerably lower than in frontal motor cortical areas controlling forelimb movements, and this rhythm was not even detectable following D2 stimulation. It is interesting to consider that the strongest NBG oscillations were detected in motor regions related to forelimb movements, which are mainly engaged in fast hyperkinetic forms of dyskinetic behaviour in this animal model [36,62]. After all treatments, NBG oscillations were markedly less pronounced in motor regions controlling trunk movements, which are engaged in slow twisting movements [36,62]. Additionally, it is worth noting that treatment with the D2 receptor agonist did not induce NBG oscillations in the trunk motor area.
Another noteworthy difference between D1 and D2 receptor stimulation is that no NBG activity was detected in the GPe after D2 stimulation, in contrast to the high detection rate equally induced by the administration of the D1 receptor agonist or L-DOPA in this region. Treatment with L-DOPA has previously been reported to induce NBG oscillations in the GPe in 6-OHDA-lesioned rats [15,16], but this is the first report comparing the effects of D1 and D2 receptor agonists on GPe oscillatory activities. We demonstrate that the effects of L-DOPA on pallidal NBG oscillations are mimicked by the D1 receptor agonist in terms of both detection rate and power, with no effects of the D2 receptor agonist, suggesting that these oscillations in GPe are clearly driven by D1 stimulation. It is interesting to note the almost identical patterns induced in all parameters of NBG activity in both GPe and M1Tr, especially considering that these regions had a tight functional coupling (discussed in further detail below). In consideration of the recently revealed detailed topography of cortico-basal ganglia-thalamic networks [63], it is worth noting that our electrode implants aimed at a rather caudal-lateral part of GPe (see Table 1). This might explain the tighter functional coupling found between GPe and the more caudally located trunk motor area compared to rostral forelimb areas.

Effects of D1 and D2 Receptor Agonists on Theta Oscillatory Activity
In conjunction with motor state transitions occurring after the three pharmacological treatments, a prominent increase of theta oscillations (5-10 Hz) was detected in several basal ganglia structures of both the intact and lesioned hemisphere. In the deeper nuclei of the basal ganglia (GPe and SNr), theta oscillations were strongly increased by D2 stimulation, whereas the effects of D1 stimulation and L-DOPA on this rhythm were less pronounced. Alterations in the spectral contents of slow oscillatory activities have previously been reported in deep basal ganglia nuclei in PD patients after dopamine replacement therapy [12,13]. Recent studies in animal models of LID have confirmed the presence of theta oscillations in the basal ganglia after treatment with L-DOPA [15][16][17]. In accordance with our findings, these studies show a prominent increase of theta oscillations in both hemispheres, particularly present during periods of dyskinesia. However, it was proposed that some of the spectral changes are most likely associated with the more active behavioural state rather than to dyskinetic motor signs per se [15,16]. Nevertheless, recording from the subthalamic nucleus in PD patients, Alonso-Frech and colleagues [12] have shown a strong relationship between theta oscillations and dyskinesia. For example, when dyskinesias were asymmetrically expressed, theta oscillations were only present in the contralateral hemisphere. Moreover, in patients with diphasic dyskinesia, theta oscillations occurred in strict temporal coincidence with the involuntary movements. In line with these data, our results show significant increases in theta oscillations during the expression of dyskinesia, the strongest effects being seen in GPe and SNr of the lesioned hemisphere after treatment with sumanirole. These data suggest an important role of D2 receptor stimulation in the induction of theta oscillations.
It remains to be established what particular form of PD dyskinesias is mimicked by D2 receptor agonists inducing pallidal theta oscillations. Whereas all current animal models of LID (as treated with L-DOPA) reproduce choreiform "on" dyskinesias [28], PD patients exhibit a larger variety of dyskinetic manifestations that may depend on a differential engagement of D1 vs. D2 receptors. Accordingly, it has been suggested that "on" and biphasic dyskinesias are differentially dependent on D1 vs. D2 receptors, respectively [64]. Taking into consideration that increases in theta activity have been reported in PD patients experiencing diphasic dyskinesias [12], it is tempting to speculate that the D2 agonist-induced AIMs in our animal model mimic diphasic dyskinesia rather than the "on" pattern. Another aspect to consider is the relative representation of fast hyperkinetic vs. slow dystonic-like components in the involuntary movements. Indeed, our recent studies in L-DOPA-naïve 6-OHDA-lesioned mice show that, in the absence of D1 stimulation, D2 receptors selectively mediate dystonic forms of dyskinesia [62]. Interestingly, an association between dystonic motor manifestations and theta oscillations has been found in primary dystonia patients using pallidal recordings (reviewed in [9]). In our animal model, sumanirole-induced dyskinetic behaviours were characterised by a predominant representation of slow components (axial AIMs) relative to fast hyperkinetic components (limb and orolingual AIMs). Moreover, the analysis of treatment-induced rotational behaviour revealed slower turning speed after treatment with sumanirole compared with SKF82958. Taken together, our data suggests that the stronger increase in pallidal and nigral theta oscillations induced by sumanirole may underlie the different pattern of dyskinetic behaviours induced by D2 compared to D1 receptor agonists.

Phase Analysis of NBG Oscillations
To study the network dynamics of NBG oscillations, we examined the changes in inter-structure connectivity during peak dyskinesia that, in parallel to intrinsic network changes induced by the dopaminergic denervation and the pharmacotherapy, may facilitate the transmission of these oscillations between different brain regions [54,65,66]. Phase analysis revealed that NBG oscillations in the M1FL were leading motor cortical and striatal structures, with a stronger effect following treatment with the D1 receptor agonist, which suggests that this region might be a key driver of these oscillations. It is, however, important to keep in mind that both cortex and the basal ganglia are under strong influence of thalamic input, which potentially could be a common pacemaker of these high-frequency oscillations in several parts of the circuit (reviewed in [23]). In line with the consistent phase differences, primary motor cortical areas (M1FL and M1Tr) showed strong functional connectivity intracortically as well as to the striatum upon D1 receptor stimulation, independent of the effects on oscillatory power. Interestingly, M1Tr and RFA also showed functional connectivity to GPe, supporting the theory that high-frequency oscillations are intrinsically generated in the cortical network, from which they spread to subcortical structures (see [23]), with a more pronounced role of the D1 receptor (as covered above). Collectively, these findings suggest a fundamental role of altered D1 receptor modulation of cortico-basal ganglia connectivity in the emergence of NBG oscillations underlying LID. Fig. 8 Graphical summary representing treatment effects on corticobasal ganglia oscillations during peak dyskinesia. Effects of L-DOPA, the D1R agonist SKF82958, or the D2R agonist sumanirole on NBG (top row) or theta oscillations (bottom row) recorded in the lesioned hemisphere during peak dyskinesia severity. The detection rate within each structure is represented using the indicated colour scale. Abbre-viations: rostral forelimb area (RFA), primary motor cortex forelimb area (M1FL), primary motor cortex trunk area (M1Tr), dorsomedial striatum (DMS), dorsolateral striatum (DLS), globus pallidus pars externa (GPe), substantia nigra pars reticulata (SNr), and narrowband gamma (NBG)