Forty-four participants were recruited, comprising 21 patients with idiopathic Parkinson’s disease (UK Brain Bank clinical diagnostic criteria) and 23 age-matched healthy controls. Four patients and two controls were excluded due to poor-compliance in-scanner with the task, or technical issues during the scanning, resulting in 17 analysed patient datasets and 21 analysed control datasets. Inclusion criteria were age 50–80 years, right-handed and no previous history of neurological or psychiatric illness apart from Parkinson’s disease. The revised Addenbrooke’s Cognitive Assessment Scale (ACE-R), including the mini-mental state examination (MMSE), was administered to all participants. No patients had clinical presentations of dementia, and there was no difference in ACE-R scores between groups (Table 1). Patients were also assessed on the Unified Parkinson’s disease rating scale.³⁹ Clinical assessments and experimental tasks were performed on usual dopaminergic medication, and dopaminergic dose equivalents were calculated.⁴⁰ Experimental protocols including written informed consent conformed to the guidelines of the Declaration of Helsinki and were approved by the Cambridge Research Ethics Committee (CRE code: 07/H0307/64).

Stimuli were presented using Matlab (MathWorks, Natick, MA, USA) and Psychtoolbox in a quiet and dimly lit room. For training, stimuli were displayed on a cathode-ray tube (CRT-type) monitor at 60 cm distance from the observer. For the scan session, stimuli were projected on a screen at 130 cm distance (60 Hz refresh rate) with equivalent pixel resolution of 0.03°. On each trial, four random dot kinematograms were displayed within four circular apertures (4° diameter) positioned along a semi-circular arc (3.4° eccentricity) on a black background. Two hundred dots were displayed during each frame and spatially displaced to introduce apparent downward motion (6°/s velocity). Motion coherence was manipulated by allowing only a certain proportion of dots to move downward on each frame whilst the rest were randomly reallocated. Motion coherence level was varied across trials. The 1.5 s long coherent motion interval was preceded and followed by intervals of zero-coherence levels lasting 1 s and 0.5 s, respectively. This was to avoid large stimulus sensory-evoked potentials elicited by abrupt stimulus onset/offset, which might mask decision processes.

A practice session adopting 100% coherent stimuli allowed participants to familiarize themselves with the task. The practice session ended once participants reached 90% accuracy across all trial types. In the following psychophysical session, motion coherence was varied between trials to estimate individual motion thresholds. Eight logarithmically spaced motion coherence levels (0 0.5 0.10. . .0.9) were used (32 trials per level). Each training session comprised 16 blocks of 32 trials. Feedback was provided for correctness of responses as well as for too early or too late responses (100 ms and 2.5 s from motion coherence onset, respectively).

The psychophysical task was combined with a control match-to-sample task where occasionally (P = 0.2) after a correct choice, participants had to compare the location of a set of grey discs with the location of the previously displayed coherent stimuli. This was to ensure that participants perceived all the available options (i.e. coherent stimuli) before committing to a decision. To report a match, participants had to press any button and withhold a response otherwise. A trial was considered as correct only when both choice and matching were correct. Trials with un-matching responses were discarded and repeated within the session.

To tailor low and high levels of perceptual uncertainty to individual motion sensitivity across number of coherent stimuli, the motion-coherence dependent accuracy of each trial type (e.g. either 1 or 3 coherent stimuli) was fitted using a maximum likelihood method, with the Log-Quick function

where α is the threshold, β is the slope and x is the coherence level. To obtain the correct proportion for each trial type, Flog was scaled by

where γ is the guess rate, and λ is the lapse rate controlling the lower and upper asymptote of the psychometric function, respectively.

Individual low and high perceptual uncertainty levels for each trial type were estimated as the 75th and 90th percentiles of the psychometric function.¹⁹

Participants performed a task in which they had to report a downward coherent stimulus by pressing the corresponding button.¹⁹ The number of available coherent stimuli defined two trial types: low action uncertainty trials, where a single coherent stimulus commanded which button to press, and high action uncertainty trials, where three coherent stimuli required the participants to press any one of the three corresponding buttons (making a ‘fresh choice, regardless of what you have done in previous trials’). Equal emphasis was placed on the speed and accuracy of the responses, and guessing was discouraged. Participants were instructed to fixate on a central red mark throughout the trial. Each trial started with the presentation of the fixation mark and stimuli onset ensued after a variable interval comprised between 0.5 s and 1 s. The imaging session was preceded by one training psychophysical session scheduled on a separate day (maximum 4 days between sessions). In the scan sessions, coherence levels were fixed to the individual motion thresholds corresponding to high and low levels of perceptual uncertainty, the match-to-sample task was removed, and no feedback was provided except for too early or too late responses. Levels of perceptual and action uncertainty were randomly interspersed across trials. Each session consisted of 10 blocks (total 640 trials per participant) separated by a short rest (total duration of the session ∼60 minutes).

Following MEEG, participants underwent structural MR imaging using a 3T Siemens Tim Trio scanner with a 32-channel phased-array head coil. A T₁-weighted magnetization-prepared rapid gradient-echo image was acquired with repetition time (TR) = 2250 ms, echo time (TE) = 3.02 ms, matrix = 192 × 192, in-plane resolution of 1.25 × 1.25 mm, 144 slices of 1.25 mm thickness, inversion time = 900 ms and flip angle = 9°. Images were first aligned to a canonical average image in MNI space, before segmentation and calculation of total intracranial volume (TIV). After segmentation, a study-specific DARTEL template was created from the patient scans and the seventeen closest age-matched controls. The remaining controls were then warped to this template. The templates were aligned to the SPM standard space and the transformation applied to all individual modulated grey-matter segments together with an 8 mm FWHM Gaussian smoothing kernel.

An Elekta Neuromag Vectorview System (Helsinki, Finland) simultaneously acquired magnetic fields from 102 magnetometers and 204 paired planar gradiometers, and electrical potential from 70 EEG electrodes (70 Ag-AgCl scalp electrodes in an Easycap—GmbH, Herrsching, Germany—extended 10–10% system). Additional electrodes provided a nasal reference, a forehead ground, paired horizontal and vertical electro-oculography, electrocardiography and right arm electromyography. All data were recorded and digitized continuously at a sample rate of 1 kHz and high-pass filtered above 0.01 Hz.

Before scanning, head shape, the locations of five evenly distributed head position indicator coils, EEG electrodes location and the position of three anatomical fiducial points (nasion and left and right pre-auricular) were recorded using a 3D digitizer (Fastrak Polhemus Inc., Colchester, VA). The initial impedance of all EEG electrodes was optimized to below 10 kΩ, and if this could not be achieved in a particular channel, or if it appeared noisy to visual inspection, it was excluded from further analysis. The 3D position of the head position indicators relative to the MEG sensors was monitored throughout the scan.

These data were used by Neuromag Maxfilter 2.2 software,⁴¹ to perform environmental noise suppression, motion compensation and Signal Source Separation. Subsequent analyses were performed using in-house Matlab (MathWorks) code, SPM12 (www.fil.ion.ucl.ac.uk/spm) and EEGLab⁴² (Swartz Center for Computational Neuroscience, University of California San Diego). Artefact rejection was performed through separate independent component analysis decomposition for the three sensor types. For EEG data, components temporally and spatially correlated to eye movements, blinks and cardiac activity were automatically identified with EEGLab’s toolbox ADJUST.⁴³ For MEG data, components were automatically identified that were both significantly temporally correlated with electro-oculography and electrocardiography data, and spatially correlated with separately acquired topographies for ocular and cardiac artefacts. A further independent component analysis identified that MEG and EEG components significantly temporally correlated with tremor-related surface electromyographic activity (2–8 Hz) recorded in Parkinson’s disease patients. Artefactual components were finally projected out of the dataset with a translation matrix.

The continuous artefact-corrected data were low-pass filtered (cut-off = 100 Hz, Butterworth, fourth order), notch filtered between 48 and 52 Hz to remove main power supply artefacts, down-sampled to 250 Hz and epoched from −1500–2500 ms relative to motion coherence onset. EEG data were referenced to the average over electrodes. MEG and EEG data were combined before inversion into source space⁴⁴ using the Minimum Norm algorithm as implemented by SPM12. Notably, combined MEG and EEG allows a better localization of neural sources than each technique on its own.⁴⁴ The forward model was estimated from each participant’s anatomical T₁-weighted MRI image. All conditions were included in the inversion to ensure an unbiased linear mapping. The source images were spatially smoothed using an 8 mm FWHM Gaussian kernel.

To decompose behavioural performance into latent variables underpinning decisions, we fitted LBA models to each participant’s reaction time data. The LBA model belongs to the broad class of accumulation-to-threshold models of decision-making but is more tractable than drift-diffusion models for n-way decisions while still remaining physiologically informative.⁴⁵ We followed the same logic of previous work¹⁹ and opted for a ‘unitary’ model where both perceptual and action uncertainty concur in determining a participant’s performance in a given trial. The unitary model accommodated the empirical data under the assumption of a constant non-decision time across experimental conditions, which provides theoretical evidence against the need of an extra stage of processing.⁴⁶ Accordingly, each accumulator linearly integrates the decision-evidence (or the intention) over time in favour of one action, and the decision is made when the accumulated activity reaches a threshold (see Fig. 1B). In our task, possible actions correspond to a button press from one of four fingers, each modelled by independent accumulators. When three valid actions are available, three accumulators are engaged with activation starting at levels independently drawn from a uniform distribution, and increasing linearly over time with an accumulation rate v drawn from an independent normal distribution. A response is triggered once one accumulator wins the ‘race’ and reaches a decision bound b. When only one action is available, only the accumulator corresponding to the available action is engaged. Predicted reaction time is given by the duration of the accumulation process for the winning accumulator, plus a constant non-decision time t₀ representing the latency associated with stimulus encoding and motor response initiation.

The LBA model was fitted to the behavioural data using a bounded version of MATLAB’s fminsearch function. We employed the Maximum Likelihood Estimation (MLE) procedure through a custom MATLAB module. To enhance numerical convergence, we adopted multi-start optimization. This method was uniformly applied to every model variant and participant dataset, with 25 optimization runs for each. Each optimization run was preceded by 100 random searches to determine the best starting parameter values. The reliability of our findings was verified by repeated optimization sessions converging on similar solutions, allowing minor numerical variations.

To improve computational efficiency and reduce multiple comparisons, we reduced the dimensionality of the MEEG data by parcellating the cortical surface into a set of 96 regions of interest (ROIs) defined using the Harvard-Oxford cortical atlas (FSL, FMRIB, Oxford). The dynamic of each ROI was represented by a single time-course, obtained by extracting the principal component from the vertices lying within the given ROI.¹⁹ The reconstructed sources within each ROI were first bandpass-filtered in either beta (13–30 Hz) or gamma (31–90 Hz) frequency bands. The coefficients of the principal component accounting for the majority of the variance of the vertices within each ROI were then taken as an appropriate representation of source activity for that region. Next, to estimate the power oscillations on a single-trial basis, we extracted frequency-specific signal envelope modulations using a Hilbert transform of the source data from each reconstructed ROI (epochs from 500 ms before to 1500 ms after coherence onset). The Hilbert’s envelope is a convenient measure of how the power of the signal varies over time in the frequency range of interest, and thus particularly suited to capture relatively slow fluctuations associated to the instantaneous accumulation of evidence/intentions. The power estimates of individual participants were down-sampled to 100 Hz and normalized by their baseline (from 400–100 ms before coherence onset).

For all the frequentist statistical analyses, the threshold of statistical significance was set to α = 0.05, corrected for multiple comparisons where appropriate.

Demographic and clinical characteristics of the participants are shown in Table 1. The groups were matched for gender (Fisher’s exact test, P∼1.00, odds ratio = 0.88), for age (Wilcoxon rank-sum test W = 120.50, P = 0.13, two-tailed), and performed similarly on cognitive tests (ACE-R Wilcoxon rank-sum test W = 161.50, P = 0.79, two-tailed).

To confirm the integrity of the network of sensory, motor and association cortices upon which our experiment relies, we used voxel-based morphometry to compare grey matter volume across the groups (Fig. 2; see Table 1 for participant characteristics). We adopted classical whole brain statistical parametric mapping t-test and Bayesian null tests to assess the presence or absence of structural differences between groups. Consistent with the early stage of the disease and lack of cognitive impairment, there was evidence for lack of atrophy in most of the cortex, and no significant evidence for atrophy anywhere in cortex.

Individual motion coherence thresholds estimated from the training session (Fig. 3) did not significantly differ between the two groups (F(1,36) = 1.59, P = 0.216, η^G2=0.037, BF = 0.92). Trials from the experimental session where responses were shorter than 100 ms or longer than 2100 ms were omitted from the behavioural and modelling analyses (controls: 15.34%, Parkinson’s disease: 20.42%). Analysis of the removed trials showed that patients did not significantly differ from controls in the number of premature responses (i.e. RT < 100 ms; Wilcoxon test: W = 397, P = 0.35, two-tailed) or omitted responses (W = 345, P = 0.06, two-tailed).

For the remaining trials, repeated-measures ANOVA showed that reaction times across groups were significantly faster for the low perceptual uncertainty condition compared with the high perceptual uncertainty condition (F(1,36) = 130.25, P < 0.001, η^G2=0.197, BF > 100) confirming the efficacy of the estimated motion coherence thresholds. We confirmed the expected lack of a significant difference between action uncertainty levels in a n-way race⁵¹ (F(1,36) = 1.99, P = 0.167, η^G2=0.005, BF = 0.648), and lack of interaction between action and perceptual uncertainty. There was no significant difference in reaction times across groups (F(1,36) = 2.46, P < 0.126, η^G2=0.053, BF = 1.08). Overall, Parkinson’s disease patients’ choices were less accurate than healthy controls (accuracy F(1,36) = 6.77, P = 0.013, η^G2=0.130, BF = 4.68); For both groups, accuracy decreased with high perceptual uncertainty (F(1,36) = 18.01, P < 0.001, η^G2=0.035, BF > 100).

Changes in the accumulation rate alone (model two, exceedance probability = 1; Fig. 1B) best accounted for the effects of uncertainty on performance in both groups (P(H=∣Data) = 0.997, BF > 100). The goodness of fit of model two (henceforth, the LBA model) was confirmed by posterior predictive checks and parameter recovery.

A mixed-design repeated-measures ANOVA on the LBA model parameters revealed slower evidence accumulation under high uncertainty levels for both action (accumulation rate: F(1,36) = 190.97, MSE=2.33, P<0.001, η^G2=0.453, BF > 100) and perceptual (accumulation rate: F(1,36) = 109.22, MSE=0.77, P<0.001, η^G2=0.136, BF > 100) manipulations. Overall, evidence accumulation was slower in the Parkinson’s disease group (F(1,36) = 4.69, MSE=11.41, P=0.037, η^G2=0.091, BF = 2.43), with a significant interaction between group and action uncertainty (accumulation rate: F(1,36) = 5.65, MSE=2.33, P=0.023, η^G2=0.024, BF = 26.95): controls showed larger changes in accumulation rates than patients in response to action uncertainty. Post hoc comparisons using Bonferroni-corrected t-tests showed that the interaction was caused by accumulation rates significantly slower in the Parkinson’s disease group compared to controls under low action uncertainty (Low action uncertainty: ΔM=1.785, t(50.10)=2.952, P=0.0048; High action uncertainty: ΔM=0.602, t(50.10)=0.995, P=0.3245).

Taken together, our results show that, regardless of group, the experimental manipulation of perceptual and action uncertainty modulated accumulation rates in line with previous reports.^2,19,50 However, compared to healthy controls, the Parkinson’s disease group was characterized by a reduced reactivity of accumulation rates to changing uncertainty levels in the action domain. Specifically, the intention to move was accumulated significantly slower than controls under low levels of action uncertainty (i.e. when a single specific action was to be taken).

The temporal evolution of the combined MEG and EEG (MEEG) power envelope from each ROI of the parcellated surface served as the signal for our analysis in beta (13–30 Hz) and gamma (31–90 Hz) bands. Confirming results from young healthy adults,¹⁹ the LBA model predictions were inversely correlated with the MEEG oscillations in beta and gamma bands (Figs 4 and 5).

Specifically, after coherence onset, neural activity desynchronized in a graded fashion and peaked approximately before response suggesting a form of threshold mechanism.^17,52 For perceptual decisions, the LBA model predicts that the accumulated decision-evidence will ramp quickly with low perceptual uncertainty, and slowly with high perceptual uncertainty.

Accordingly, desynchronization of beta power-envelopes averaged across trials and ROIs was larger (P = 0.0004 Bonferroni-corrected, cluster-based permutation test) for low than high perceptual uncertainty^19,52 in controls (Fig. 4). Such a trend was seen in the Parkinson’s disease group (P = 0.088 Bonferroni-corrected). When a response is chosen between multiple options, the race underlying the selection of each alternative is characterized by a larger amount of decision-evidence summed across all the racing accumulators by the time of response.⁵¹ Accordingly, desynchronization of beta power-envelopes averaged across trials and ROIs was larger for high than low action uncertainty in both groups (P = 0.0004, Bonferroni-corrected, cluster-based permutation test). Gamma power-envelopes showed a similar trend in the control group, but the effects were statistically insignificant (Fig. 5). Decision-related dynamics expressed by beta desynchronization were distributed across a wide network (Fig. 6A, mean sign-test Z across significant ROIs, controls: Z = −5.1 ± 0.25, P < 0.0001; patients: Z = −3.99 ± 0.58, P < 0.0001; FDR-corrected) similar to previous reports.^9,19

Comparisons between Z-transformed correlation values from each of the four levels of our manipulations in isolation confirmed that the quality of fit and results did not vary across groups and trial types (P > 0.05, FDR-corrected). In the gamma band, we observed a more localized mosaic of ROIs. In the control group, significant ROIs included contralateral motion sensitive areas (inferior lateral occipital region), bilateral extrastriate areas and bilateral frontal regions (comprising frontal pole and superior middle gyrus), ipsilateral motor and supplementary motor area; mean across significant ROIs: sign-test Z = −3.22 ± 0.58, P = 0.0031, FDR-corrected). In comparison, fewer ROIs survived statistical test in Parkinson’s disease patients (sign-test Z = −2.62 ± 0.31, P = 0.0084, FDR-corrected) and none of them in the left dorsal path, which is of primary interest for this study. Therefore, from now onwards, the analyses will focus only on beta frequencies.

To explore the moderation effects of disease on beta desynchronization, we fitted a linear mixed effect model to predict power amplitude with Group (controls, patients), PU (low, high) and AU (low, high) as fixed effects, and specified nested subjects and ROIs as random factors. The model’s total explanatory power was substantial (conditional R2 = 0.82). The analysis confirmed the effects of perceptual (β = 1.30, confidence interval [CI] [1.14, 1.46], P < 0.001) and action (β = −2.73, CI [−2.90, −2.57], P < 0.001) uncertainty. A significant interaction between Group and AU indicated lower reactivity of beta power to varying levels of action uncertainty in the Parkinson’s disease group than in controls (Fig. 7A; β = 0.42, CI [0.18, 0.67], P < 0.001, post hoc: controls High-Low = −2.73, SE = 0.083, t-ratio = −33, P < 0.0001; patients High-Low = −2.31, SE = 0.092, t-ratio = −25, P < 0.0001; Bonferroni-corrected). Finally, a non-parametric Wilcoxon rank-sum test confirmed exaggerated beta power (i.e. reduced desynchronization) in Parkinson’s disease compared to controls (dorsal path—both hemispheres: W = 119, P = 0.0416, left tailed: controls < patients).

To examine the decision-evidence accumulation over space and time, the onset of evidence accumulation across ROIs was identified by optimizing the split of the non-decision time before and after the accumulation period using Spearman correlation to the single-trial MEEG power envelope (Fig. 6A, see ‘Materials and methods’ for details). By tracing the spectrally resolved temporal evolution of decision onset through the visuomotor hierarchy, we found that decision-evidence accumulation emerges with distinct spatiotemporal profiles between healthy controls and Parkinson’s disease patients (Fig. 6A–C). In the control group, we replicated the results from Tomassini et al.¹⁹ Specifically, we show that in the contralateral (i.e. left) hemisphere, the beta-mediated evidence accumulation unfolds in a bottom-up cascade proceeding from caudal sensorial regions to rostral executive areas (Fig. 6A and C).

In people with Parkinson’s disease, accumulation in the beta frequency begins ‘earlier’ than in controls (∼240 ms from coherence onset; Wilcoxon rank-sum: W = 7977, P < 0.001). Crucially, the caudo-rostral beta gradient of evidence accumulation is abolished, with frontal regions beginning to accumulate evidence in parallel with caudal visual areas.

The difference between groups is shown in Fig. 6. We fitted a regression model to the mean latencies of ROIs located along the dorsal path for visuomotor decisions^53,54 (Fig. 6C). For controls (Fig. 6C left top-bottom panels), there is a gradient from caudal to rostral regions (R² = 0.44, P = 0.008), which contrasts the lack of gradient in Parkinson’s disease patients (R² = −0.08, P = 0.733). This indicates early decision processes begin in the frontal pole and the middle frontal gyrus, preceding (∼150 ms) onsets in the occipital pole.