Whole-Head Functional Near-Infrared Spectroscopy as an Ecological Monitoring Tool for Assessing Cortical Activity in Parkinson’s Disease Patients at Different Stages

Functional near-infrared spectroscopy (fNIRS) is increasingly employed as an ecological neuroimaging technique in assessing age-related chronic neurological disorders, such as Parkinson’s disease (PD), mainly providing a cross-sectional characterization of clinical phenotypes in ecological settings. Current fNIRS studies in PD have investigated the effects of motor and non-motor impairment on cortical activity during gait and postural stability tasks, but no study has employed fNIRS as an ecological neuroimaging tool to assess PD at different stages. Therefore, in this work, we sought to investigate the cortical activity of PD patients during a motor grasping task and its relationship with both the staging of the pathology and its clinical variables. This study considered 39 PD patients (age 69.0 ± 7.64, 38 right-handed), subdivided into two groups at different stages by the Hoehn and Yahr (HY) scale: early PD (ePD; N = 13, HY = [1; 1.5]) and moderate PD (mPD; N = 26, HY = [2; 2.5; 3]). We employed a whole-head fNIRS system with 102 measurement channels to monitor brain activity. Group-level activation maps and region of interest (ROI) analysis were computed for ePD, mPD, and ePD vs. mPD contrasts. A ROI-based correlation analysis was also performed with respect to contrasted subject-level fNIRS data, focusing on age, a Cognitive Reserve Index questionnaire (CRIQ), disease duration, the Unified Parkinson’s Disease Rating Scale (UPDRS), and performances in the Stroop Color and Word (SCW) test. We observed group differences in age, disease duration, and the UPDRS, while no significant differences were found for CRIQ or SCW scores. Group-level activation maps revealed that the ePD group presented higher activation in motor and occipital areas than the mPD group, while the inverse trend was found in frontal areas. Significant correlations with CRIQ, disease duration, the UPDRS, and the SCW were mostly found in non-motor areas. The results are in line with current fNIRS and functional and anatomical MRI scientific literature suggesting that non-motor areas—primarily the prefrontal cortex area—provide a compensation mechanism for PD motor impairment. fNIRS may serve as a viable support for the longitudinal assessment of therapeutic and rehabilitation procedures, and define new prodromal, low-cost, and ecological biomarkers of disease progression.


Introduction
In recent decades, functional near-infrared spectroscopy (fNIRS) has received increasing interest for non-invasively assessing cortical hemodynamic activity [1]. This optical imaging technique allows the direct measurement of either variations in (i.e., ∆[HbO 2 ] and ∆[HbR]) or absolute concentrations of oxygenated HbO 2 and deoxygenated (i.e., reduced)

Group-Level Activation Maps
Group-level activation maps across chromophores and task conditions for the ePD and mPD groups are separately represented in Figures 1 and 2, while ROI-based significant activations are reported in Table 2. In Section 4.6, we report the list of selected ROIs and their correspondence with functional and anatomical cortical areas. (Top) Significant group-level activation maps referring to right-grasp task. Results were corrected for multiple comparisons (pFDR < 0.05). (Bottom) ROI analysis of group-level right-grasp task for the ePD and mPD groups: (*) significant contrasts not corrected for multiple comparisons (p < 0.05); (**) significant contrasts corrected for multiple comparisons (pFDR < 0.05). The correspondence between ROIs (i.e., horizontal axis) with functional and anatomical cortical areas is reported in Section 4.6. (Top) Significant group-level activation maps referring to right-grasp task. Results were corrected for multiple comparisons (p FDR < 0.05). (Bottom) ROI analysis of group-level right-grasp task for the ePD and mPD groups: (*) significant contrasts not corrected for multiple comparisons (p < 0.05); (**) significant contrasts corrected for multiple comparisons (p FDR < 0.05). The correspondence between ROIs (i.e., horizontal axis) with functional and anatomical cortical areas is reported in Section 4.6.

Figure 2.
(Top) Significant group-level activation maps referring to left-grasp task. Results were corrected for multiple comparisons ( < 0.05). (Bottom) ROI analysis of group-level left-grasp task for the ePD and mPD groups: (*) significant contrasts not corrected for multiple comparisons ( < 0.05); (**) significant contrasts corrected for multiple comparisons ( < 0.05). The correspondence between ROIs (i.e., horizontal axis) with functional and anatomical cortical areas is reported in Section 4.6. (Top) Significant group-level activation maps referring to left-grasp task. Results were corrected for multiple comparisons (p FDR < 0.05). (Bottom) ROI analysis of group-level left-grasp task for the ePD and mPD groups: (*) significant contrasts not corrected for multiple comparisons (p < 0.05); (**) significant contrasts corrected for multiple comparisons (p FDR < 0.05). The correspondence between ROIs (i.e., horizontal axis) with functional and anatomical cortical areas is reported in Section 4.6.
From these results, we concluded, as expected for both ∆[HbO 2 ] and ∆[HbR], that contralateral activations were the most significant ones; i.e., the left (L-SMN) and right (R-SMN) motor areas due to right-grasp (RG) and left-grasp (LG) tasks, respectively. A significant ipsilateral deactivation of the L-SMN due to LG was also found for the mPD group (t = −2.784, p = 0.007). In addition, we found that L-VIS2 for ∆[HbO 2 ]-RG, L-PFC2 for ∆[HbO 2 ]-LG and ∆[HbR]-RG, L-PFC1 for ∆[HbO 2 ]-LG, ∆[HbR]-LG, and ∆[HbR]-RG were major non-motor regions that presented significant ancillary activation. Isolated significant activations were found in the ePD group in R-VIS1 for ∆[HbR]-RG (t = −2.046, p = 0.043), L-PFC2 (t = −2.510, p = 0.013) and R-PFC2 (t = −4.678, p < 0.001) for ∆[HbR]-LG, R-PFC2 for ∆[HbR]-RG (t = −5.069, p < 0.001), R-PFC2 (t = 2.653, p = 0.009), and R-PFC1 (t = 2.253, p = 0.026) for ∆[HbO 2 ]-RG. In contrast, significant isolated activations were found in the mPD group in R-PFC2 for ∆[HbO 2 ]-LG (t = 2.024, p = 0.045), L-VIS1 (t = −3.348, p = 0.001), and L-VIS2 (t = −2.518, p = 0.013) for ∆[HbR]-LG. The ePD and mPD activation maps have been separately represented in Figures 1 and 2, and significant differences were analyzed, as shown in Table 1. Next, the mapping of ePD vs. mPD contrasts were analyzed. Obviously, the overall mapping similarities and the emphasis of motor areas were canceled, as only the contrasts were represented. Nevertheless, the patterns of significant differences were outlined and discussed. The contrast maps were extracted by a two-tailed t-test contrast, regarding task condition vs. group interactions ( Figure 3 and Table 3). LG; t = −3.706, p < 0.001 for Δ[HbR]-RG). Further details regarding significant group-level contrasts among the ePD, mPD, and ePD vs. mPD subdividisons c across ROIs, task conditions, and chromophores are provided in supplementary data (sheet groupAnalysisEPD, groupAnalysisMPD, and groupAnalysisEPDvsMPD, respectively).   We found a unilateral hemispheric difference in both ∆[HbO 2 ] (L-SMN: t = 2.228, p = 0.027) and ∆[HbR] activity (L-SMN: t = −2.120, p = 0.036) due to RG, meaning that the ePD group presented a greater contralateral activation in RG compared to the mPD group. Conversely, LG caused a significant ∆[HbO 2 ] difference in the ipsilateral left motor areas (L-SMN: t = 2.768, p = 0.006), which was counter-intuitive, as we expected to find this difference in the contralateral hemisphere was due to the hand performing the task. However, Figure 3 suggests a difference in the contralateral hemisphere due to LG, and ROI analysis revealed that there was only a slight tendency toward significance (R-SMN: t = 1.516, p = 0.132), as a comparable number of positive vs. negative t-value channels were averaged simultaneously.
Concerning the contrast between the ePD group and the mPD group in non-motor areas, most of differences were found for the LG condition and the ∆ ROIs, task conditions, and chromophores are provided in supplementary data (sheet groupAnalysisEPD, groupAnalysisMPD, and groupAnalysisEPDvsMPD, respectively).

ROI-Based Correlation Analysis
The last analysis involved the ROI-based correlation analysis (ROI-CA) between subject-level fNIRS data (i.e., beta-GLM coefficients) and clinical variables. Table 4 provides a summary of all combinations of ROI, grasping tasks (LG or RG), and chromophores (∆[HbR] or ∆[HbO 2 ]) that presented at least one significant correlation with any clinical variable. No significant difference was found for LG and ∆[HbO 2 ] combinations and, accordingly, no line is reported in Table 4. Further details regarding the full tables of correlations and partial correlations are provided in supplementary data (sheet correlationAnalysis and partialCorrelationAnalysis). Table 4. Significant results of ROI-based correlation analysis between ROI averaged data and patient characteristics at subject level (p < 0.05). Results are provided by means of Spearman's correlation coefficient ρ S (p value ). Bold highlights the significant comparisons. It is worth noting that age across ROIs and tasks was always positively correlated with ∆[HbR] in motor areas that were contralateral to the hand performing the task, positively correlated with the frontal cortex, and negatively correlated with the visual cortex. As noted above, a positive ρ S -coefficient for ∆[HbR] provides evidence of lower activation with age. Congruently, a negative ρ S for ∆[HbO 2 ] was found, although it was not statistically Given the age dependency and the unbalanced ages in the ePD and mPD groups, we decided to investigate other variables according to partial correlation, using age as the covariate and noting which variables remained statistically significant (Table 5) Table 5. Significant results of ROI-based partial correlation analysis between ROI averaged data and patient characteristics at subject level (p < 0.05), using age as covariate. Results are provided by means of Spearman's correlation coefficient ρ S (p value ). (*) Correlation values were statistically significant when not employing age as covariate. Bold highlights the significant comparisons. L-PFC1 (ρ S = 0.354, p = 0.034; no significant partial correlation) for ∆[HbR]-RG, while with R-PFC2 (ρ S = 0.394, p = 0.017; ρ Spart = 0.381, p = 0.024) and R-PFC1 (ρ S = 0.398, p = 0.016; ρ Spart = 0.399, p = 0.018) for ∆[HbR]-LG. Again, it is worth noting that the SCWT scores were significantly associated with contralateral areas to the hand performing the task (i.e., L-PFC2 and L-PFC1 with ∆[HbR]-RG; R-PFC2 and R-PFC1 with ∆[HbR]-LG), while SCWE to the ipsilateral hemisphere (i.e., R-VIS1, R-PFC2 and R-PFC1 with ∆[HbR]-RG).

fNIRS Results' Interpretation
In this study, we used fNIRS for an ecological assessment of cortical activation in PD patients at different stages. Patients performed the same motor grasping task, unveiling different activation patterns over both motor and non-motor networks. Even if their illness condition spanned from early to moderate, we found a significant difference between the ePD group and the mPD group and an association between the level of cortical activation with the clinical variables. The overall results suggest that PD motor impairment affects the activation level of motor areas, depending on the stage of the pathology, while the recruitment of non-motor areas follows different trends of activation.
Specifically, we first verified that (despite there being no significant differences in cognitive reserve and executive functions, as measured by the SCW test) the mPD group is associated with higher levels of motor impairment and disease progression, compared with the ePD group (Table 1). Then, considering the comparison of the ePD and mPD group-level activation maps, we found that the ePD group had greater levels of activation and inhibition across task conditions and ROIs. Specifically, the ePD group presented greater activation levels in the SMN and VIS areas, while the mPD group required increased activation of the PFC areas ( Figure 3 and Table 3).
This result is also supported by the last ROI-CA analysis performed for all patients, which reported significant correlations of fNIRS data with the CRIQ, disease duration, the UPDRS, and the SCW performances in the PFC and VIS areas (Table 4). Indeed, these results indicate that a worsening of clinical conditions is associated with a lower activation and inhibition of these areas (i.e., positive correlations are mostly found for ∆[HbR], while negative correlations are mostly found for ∆[HbO 2 ]). This result is corroborated by partial correlations, separating the effect of age as the covariate, which showed almost no change in their level of statistical significance (Table 5), supporting a specific effect of disease severity.

Discussion
The ePD group presented greater activation patterns in SMN areas than the mPD group, according to both LG and RG tasks. Considering that both groups performed the motor grasping task correctly, the lower levels of activation in the mPD group may provide an overall index of loss of function. Indeed, it is known that cerebral perfusion and cerebrovascular reactivity of motor regions in PD are significantly correlated with disease severity [40]. As a consequence, the motor region in PD patients tends to hypoactivation, according to disease progression and the stage of PD progression over time. Interestingly, we noticed a unilateral hemispheric difference in L-SMN due to RG (i.e., significant differences were found only in the left hemisphere), while LG presented significant differences on bilateral SMN (Figure 3 and Table 3). Nevertheless, significant differences due to RG involved additional channels in L-SMN, compared with LG. This effect could possibly be due to the right-handedness of almost all patients (Table 1), causing a greater activation of L-SMN because of the dominant hand. In addition, no significant correlations were found for the LG condition for ∆[HbO 2 ] across all clinical variables (Tables 4 and 5).
Conversely, a higher PFC activation in the mPD group, compared with the ePD group, is in line with current fNIRS studies in PD, which have investigated PFC activation during walk conditions. The shared hypothesis is that PFC compensates for motor impairment and acts as an overall index of cognitive load. Among these works, Maidan et al. [16] indicated that PFC activation levels can distinguish between different typologies of FoG, as one of the most common disturbances in PD. The role of PFC as a motor-cognitive compensation mechanism was also confirmed in a successive work considering PD without FoG [17]. In addition, a recent study concluded that PD patients with FoG present a greater PFC activation, compared with healthy controls (HC) and PD without FoG, during a lower-limb motor task, and a positive correlation with the severity of FoG [15]. Usual and dual-task (DT) walking in patients with and without FoG were also assessed by Dagan et al. [18] and Orcioli-Silva et al. [19], respectively. Patients were evaluated before and after anti-Parkinsonian medication and the results provided additional insights into the interplay between PFC activation and gait parameters. Other fNIRS studies found an increase of PFC activation during walking, DT walking, and postural stability tasks with respect to HC [20,21,29]. Cornejo et al. [24] reported that walking at an externally paced rhythm can further promote PFC activity with respect to usual walking. Similarly, Maidan et al. [26] suggested that DT walking, in contrast with obstacle negotiation, causes different involvements of PFC in PD and HC. The same authors also performed a longitudinal randomized control trial that demonstrated that a combined virtual reality and a treadmill training program helped to reduce PFC activation during obstacle negotiation and DT walking [25]. Another study, by Hoang et al. [28], revealed reduced dorsolateral-PFC activity and better gait parameters after 5 weeks of an exercise-based training program.
Moreover, our results of increased PFC activation in the mPD group, compared with the ePD group, paired to a reduction in activation in the VIS and SMN areas, may be attributed to the loss of function of these areas associated with disease progression. Indeed, recent studies indicated that early PD patients, compared with HC, presented a reduction in functional connectivity and cerebral blood flow in the SMN, VIS1, and VIS2 areas. Therefore, it is proposed as prodromal biomarker of neurovascular coupling impairment and disease progression [41]. Moreover, other evidence from MRI biomarkers indicated that motor rehabilitation in PD promotes the restoration of movement automaticity and the involvement of a frontal-parietal compensatory mechanism [42].
Among the limitations of this study, we identified an unbalanced sample size in the ePD and MPD groups, which requires caution in the generalization of the results to a broader population. Moreover, the significant age difference between the groups is another potential issue, as it is known that normal aging affects the level of automaticity during the execution of sequential movements [43]. Nevertheless, in our ROI-CA, we noticed that a partial correlation between the cortical activation level of all patients with the clinical variables, using age primarily as a covariate, did not affect the significance level. Hence, in our study, we concluded that the disease stage represented the primary driving reason for cortical activation differences in the ePD and MPD groups.

Participants
This study considered the analysis of 39 PD patients (age 69.0 ± 7.64, 38 right-handed) who were enrolled in the baseline assessment of the SIDERAˆB project [44]. The inclusion criteria were as follows: • age between 18 and 85 years (adult and older adult); • agreement to participate, with signature on the informed consent form; • clinical diagnosis of PD, according to the Movement Disorder Society (MDS) criteria [45] and disease staging between 1.5 and 3 on the Hoehn and Yahr scale [46] The exclusion criteria were as follows: • presence of comorbidities that might determine clinical instability (i.e., severe orthopedic or severe cognitive deficits); • overlapping between PD and other neurological pathologies or by PD with severe psychiatric complications, based on a pathological score in a screening test for cognitive impairment (Montreal Cognitive Assessment test-MoCA test < 17.54 [47]) Patients were subdivided into two different groups according to the different stage of the pathology, as assessed by the Hoehn and Yahr (HY) scale. Specifically, patients with HY scores of 1 or 1.5] were assigned to the ePD group, while patients were assigned to the mPD group if they had HY scores of 2, 2.5, or 3]. All patients provided their informed consent, and the study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of IRCCS Fondazione Don Carlo Gnocchi.

Clinical Charaterization and Neuropsychological Assesment
The overall patient demographics and clinical variables are reported in Table 6, while further details of the ePD group and. the mPD group are presented in the Results Section in Table 1. Among the clinical variables, we employed the MDS-UPDRS and disease duration to determine if the ePD and mPD groups significantly differed in terms of motor impairment, while the CRIQ [48] was used to verify whether the subgroups differed in terms of cognitive reserve and, hence, whether measured hemodynamic variations by fNIRS could be entirely attributed to PD pathology. Further variables, such as EQ5D5L [49] and education, were also employed for participants' characterization. Finally, a computerized version of the SCW test by Caffarra et al. [50] was administered as a neurophysichological test score to evaluate the activation of the subjects' non-motor areas. This test was administered on the same day and in correspondence with the fNIRS acquisition (as described in Section 4.3). All patients were in the ON-state of the anti-Parkinsonian medications when the SCW test was performed.
This test was subdivided into three trials, indicated as color (C), word (W), and interference trials (CW). The former two trials required patients to read the name of color words (i.e, "red", "blue", and "green") and to state the color of squared shapes. Conversely, in the CW trial, color words were dyed with a hue that did not represent the semantic meaning of the word, thus requiring the patients to state the hue instead of the color word (i.e., if the color word "red" was presented in blue ink, the participant had to respond "blue"). The computerized version required patients to respond to these trials by pressing three color-coded buttons. The resulting SCWE and SCWT scores were computed, respectively, as the difference between the number of errors or the time of the CW trial and the average of the C trial and the W trial.
All demographic data and clinical variables were tested for normality by the Shapiro-Wilk test for the whole group and for the ePD and mPD subdivisions in Jamovi 1.8.1.0 (the Jamovi project (2021), https://www.jamovi.org), allowing us to define appropriate statistical tests for the group's characterization (Section 2.1) and to carry out ROI-based correlation analysis (Sections 2.3 and 4.6). We also tested for significant differences in ePD and mPD group demographics and clinical variables, according to either an independent sample t-test or a Mann-Whitney U-test, depending on the data distribution (see Table 1 in the Results Section).

fNIRS Assessment
The experimental data were acquired from IRCCS Fondazione Don Carlo Gnocchi, Milan. We employed the NIRScoutX 32 × 32 (NIRx Medizintechnik, Berlin, Germany) to provide functional measurements at 760 and 850 nm. This continuous wave (CW) system used 32 LEDs sources and 32 avalanche photodiode detectors, which resulted in 102 measurement channels (i.e., source-detector combination pairs). The resulting sampling frequency was 1.9531 Hz. Sources and detectors were placed over the scalp positions of each subject, according to the International 10/5 system [51], thus identifying the channel positions as the midpoint between the source and detector positions (i.e., in the supplementary data at optodeConfigurationROI sheet, we state the correspondence between the source and detector number with International 10/5 locations).
In line with neuropsychological testing by SCW, all patients were in the ON state of the anti-Parkinsonian medications during the fNIRS acquisition. Because PD is a neurodegenerative disorder that is mainly characterized by motor deficits, a block-design motor-grasping task was employed to elicit the activation of the motor areas [52] and to observe ancillary activity in the non-motor areas. Participants were asked to sit still in a dimly lit room and to squeeze repeatedly two sponge balls that were placed within the palms of their left and right hands. A stimulus presentation PC running E-Prime 3.0 (Psychology Software Tools, Pittsburgh, PA, USA) was synchronized with the fNIRS system and placed in front of each participant. The resting condition was represented by a fixation cross placed in the center of the screen. The task condition involved the blinking at 0.8 Hz of the cross on either the left or right side of the screen, requiring participants to follow this blinking by pressing and releasing the ball with left or right hand. The stimuli presentation was repeated 10 times per side and randomized in left vs. right order. The task condition lasted 10 s and was followed by 20 s of resting.
Two experienced fNIRS users visually examined the performance of the participants to control the correct execution of the task, observing no major differences among them. No instrumental assessment of the motor performance was employed, as the task only required the participants to press and release a simple rubber ball.

Pre-Processing and Subject-Level Statistics of the fNIRS Data
A conceptual framework of the analyses performed in this work is provided in Figure 4. The pre-processing and statistical analysis of fNIRS data is still an open field in fNIRS literature and the procedure is non-standardized due to the heterogeneity of fNIRS applications, instrumentation, and experimental paradigms [53][54][55][56]. The estimation of fNIRS hemodynamic response assumes an even more important role when considering applications contextualized in clinical domains. In a previous work, we carefully considered these aspects by comparing the most-used artifact reduction algorithms in CW-fNIRS, without employing auxiliary sensing systems, as well as quantifying the reduction in signal energy at each step of a processing pipeline [57].

Group-Level Statistics of fNIRS data
A linear mixed-effect model (LMEM) was computed for the GLSA of fNIRS data using NIRS Brain AnalyzIR toolbox. This model separately tested the effect of LG and RG conditions (i.e., brain activity) for each group, while controlling for subject as random effect. Accordingly, this model was expressed according to Wilkinson-Rogers notation as ~ 1 : (1| ), where indicated the GLM regression coefficients deriving from SLSA, referred to group membership (i.e., ePD vs. mPD), referred to the subject identifier, (1| ) was a random intercept for each subject, and : was the interaction between condition and group terms. The results of the LMEM allowed us to compute three different typologies of grouplevel activation maps. First, we separately estimated the activation due to either LG or RG conditions in the ePD and mPD groups. Positive indicated activation, compared with the baseline condition (see the results in Figures 1 and 2). Then, a two-tailed t-test contrast was applied to LG-ePD interaction vs. LG-mPD interaction, and similarly to RG-ePD interaction vs. RG-mPD interaction. Positive indicated that the mean activation level in the ePD group was higher than in the mPD group (see the results in Figure  3). In all cases, channel-wise activation was considered statistically significant if the resulting t was associated with < 0.05. The false-discovery-rate method was still employed for multiple comparisons correction.

ROI-Based Correlation Analysis between fNIRS and Clinical Variables
Group-level activation maps allowed us to carry out numerical and visual evaluations of which channels were significantly activated or de-activated during LG or RG. A further refinement was provided by an ROI analysis addressing cortical areas, which required a correspondence between the whole-head channel configuration and functionalanatomical areas. In this work, we translated our previous considerations over a clinical dataset. Data were analyzed in Matlab R2018b (The MathWorks, Inc., Natick, MA, USA) through userdefined scripts integrating Homer2 [58] and NIRS Brain AnalyzIR [59] toolboxes. Before data acquisition, two experienced fNIRS users verified that all channels presented a reasonable coupling with scalp positions. This information was provided as a quality measure of the signal-to-noise ratio by the acquisition software of the employed fNIRS system.
Then, raw data were converted into optical density variations, corrected by motion artifacts through a wavelet-based method [60], and high-pass filtered with a 0.01 Hz cutoff frequency. Then, we applied principal component analysis to remove up to 80% of data variance, thereby reducing the physiological interference derived from slow-varying oscillations and scalp-related superficial confounds, such as heart rate, Mayer waves, and respiratory oscillations [61,62]. We approached the subject-level statistical analysis (SLSA) of the fNIRS data according to the general linear model (GLM). The GLM parameter estimation required us to avoid lowpass smoothing filtering (added to the mentioned high-pass filtering) to prevent artifactual temporal correlations and to obtain white noise residuals [63,64]. We specifically employed the autoregressive iteratively reweighted least squares (AR-IRLS) by Barker et al. [65]. In brief, this algorithm adapts the traditional GLM analysis by iteratively estimating the coefficients of an AR(p) model (where p is the order of the model) to apply pre-whitening and reduce serially correlated errors, which are typical in fNIRS data. Additionally, the design matrix employed a constant regressor and the convolution between boxcar functions (the period is set as stimuli duration) and a canonical hemodynamic response function (cHRF) derived from fMRI literature [66] (with a maximum peak at 6 s). Convolution by the first derivative of cHRF was applied as a further regressor to account for the temporal dispersion of the real hemodynamic response, compared twitho cHRF.
Prior to group-level statistical analysis (GLSA), subject-specific GLM data were also channel-wise contrasted for task vs. resting conditions according to a two-tailed t-test. This approach allowed us to reduce the memory computation required by the group-level linear mixed-effects model (Section 2.4), while simultaneously controlling for the statistical significance of fNIRS data. Specifically, the 10 s LG and RG conditions were compared, respectively, to the following 20 s of resting periods. For the sake of simplicity, in the group-level analysis (Section 4.5), we indicated for LG and RG the contrast between task and resting conditions.
Resulting channel-wise positive t value associated to p < 0.05 indicated that activation occurred at that given channel, while negative t value indicated inhibition. The false discovery rate by the Benjamini-Hochberg method was also used for multiple comparisons correction.

Group-Level Statistics of fNIRS data
A linear mixed-effect model (LMEM) was computed for the GLSA of fNIRS data using NIRS Brain AnalyzIR toolbox. This model separately tested the effect of LG and RG conditions (i.e., brain activity) for each group, while controlling for subject as random effect. Accordingly, this model was expressed according to Wilkinson-Rogers notation as beta ∼ −1 + group : cond + (1|subject), where beta indicated the GLM regression coefficients deriving from SLSA, group referred to group membership (i.e., ePD vs. mPD), subject referred to the subject identifier, (1|subject) was a random intercept for each subject, and group : cond was the interaction between condition and group terms.
The results of the LMEM allowed us to compute three different typologies of grouplevel activation maps. First, we separately estimated the activation due to either LG or RG conditions in the ePD and mPD groups. Positive t value indicated activation, compared with the baseline condition (see the results in Figures 1 and 2). Then, a two-tailed t-test contrast was applied to LG-ePD interaction vs. LG-mPD interaction, and similarly to RG-ePD interaction vs. RG-mPD interaction. Positive t value indicated that the mean activation level in the ePD group was higher than in the mPD group (see the results in Figure 3). In all cases, channel-wise activation was considered statistically significant if the resulting t value was associated with p < 0.05. The false-discovery-rate method was still employed for multiple comparisons correction.

ROI-Based Correlation Analysis between fNIRS and Clinical Variables
Group-level activation maps allowed us to carry out numerical and visual evaluations of which channels were significantly activated or de-activated during LG or RG. A further refinement was provided by an ROI analysis addressing cortical areas, which required a correspondence between the whole-head channel configuration and functional-anatomical areas.
Specifically, the ROI analysis consisted of averaging the statistical results at SLSA and GLSA according to predefined lists of source-detector pairs placed in correspondence with functional-anatomical areas. We considered as functional cortical parcellations the bilateral Brodmann areas referred to the Colin27 atlas [67], due to its wide adoption in fNIRS applications [68] and software [9,58,59]. Specifically, we considered the composition of multiple bilateral Brodmann areas: BA1-2-3-4 as the primary sensorimotor network areas (left L-SMN and right R-SMN), BA17 and BA18-19 as the primary and secondary visual areas, respectively, in the occipital cortex (L-VIS1 and R-VIS1 vs. L-VIS2 and R-VIS2), and BA46-9 and BA45-47 as the high-level and low-level sites of executive functions in the dorsolateral-and ventrolateral-prefrontal cortex areas (L-PFC1 and R-PFC1 vs. L-PFC2 and R-PFC2).
We employed the Nirstorm package in Brainstorm [69] to virtually register the probe configuration in the Colin27 atlas. Then, we defined each ROI-specific list of source-detector pairs by computing the sensitivity profile of the channels to cortical anatomy via Monte Carlo simulations. Channels that presented a significant contribution to the mentioned cortical parcellations (i.e., channels that collected at least 20% of the ROI signal) were included in the respective ROI list of source-detector pairs. This analysis would possibly be applied to subject-specific anatomies to properly identify cortical parcellations referred to measurement channels. However, it is reported that an atlas-based approach still leads to reasonable results in image reconstruction when anatomical MRI is unavailable [70]. A summary table of functional vs. anatomical areas correspondence is provided in Table 7, while further information regarding source-detector correspondences to ROIs is provided in the supplementary materials (sheet optodeConfigurationROI). Finally, an ROI-based correlation analysis (ROI-CA) was conducted to identify whether there was a relationship between subject-specific cortical activations-as measured by GLM beta-coefficients from SLSA-and patient characteristics. More precisely, we were interested in assessing a relationship between the motor-grasping task with either motor regions (i.e., where we expected to find the most prominent activation) or non-motor regions (i.e., where we expected to find ancillary significant brain activity).
ROI-CA involved two successive steps. First, significant channels in the GLSA results, contrasted for differences in the ePD group and the mPD group, were noted for each selected ROI. Then, an ROI average was applied to the SLSA results, contrasted for LG and RG tasks and baseline by limiting the number of channels to the ones identified in the first step. This approach provided a data-driven method to restrict the correlation analysis of the fNIRS data for brain areas, primarily contributing to activation across all subjects and reducing the possibility of averaging non-significant activations at subject level.
We applied Spearman's correlation in Jamovi 1.8.1.0, due to data distribution (see Table 1 in the Results Section), providing a ρ S -correlation coefficient for each item of fNIRS data to the clinical variable relationship. Because we were interested solely in the association of patients' characteristics with brain activity, we performed this ROI-CA for all subjects, including both the ePD group and the mPD group. Among the clinical variables listed in Table 6, we considered the correlation of fNIRS data with age, the CRIQ, disease duration, the UPDRS, the SCWE scores, and the SCWT scores. Indeed, the correlation with age, the CRIQ, disease duration, and the UPDRS allowed us to understand whether variations in the demographics and clinical outcomes were related to variations in the levels of brain activity. Conversely, the correlation of SCWE and SCWT scores provided more insights into significant ancillary activation of the non-motor areas. Finally, as age presented a significant correlation with both the motor areas and the non-motor areas, we also applied a partial Spearman's correlation using age as a covariate, thereby noticing whether previous significant correlations were preserved.

Conclusions
This work may open new scenarios in the investigation of the progression of PD motor and non-motor impairment, using fNRIS as a non-invasive and ecological neurophysiological correlate of cortical activity. fNIRS data, if used in conjunction with current clinical outcomes, will provide viable support for the monitoring of disease progression and for the longitudinal assessment of therapeutic and rehabilitation procedures. In this study, we approached the aspect of PFC compensation for motor impairment from a different standpoint by using a grasping task instead of a walking task as in the current fNIRS literature. Data from our fNIRS signals show significant differences in the frontal-occipital regions and reveal that VIS areas are involved in this mechanism. Taken together, our results demonstrate that fNIRS could provide viable support for the longitudinal assessment of therapeutic and rehabilitation procedures and define new prodromal, low-cost, and ecological biomarkers of disease progression.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data are not publicly available due to privacy issues.