Evidence of Neurovascular Un-Coupling in Mild Alzheimer’s Disease through Multimodal EEG-fNIRS and Multivariate Analysis of Resting-State Data

Alzheimer’s disease (AD) is associated with modifications in cerebral blood perfusion and autoregulation. Hence, neurovascular coupling (NC) alteration could become a biomarker of the disease. NC might be assessed in clinical settings through multimodal electroencephalography (EEG) and functional near-infrared spectroscopy (fNIRS). Multimodal EEG-fNIRS was recorded at rest in an ambulatory setting to assess NC and to evaluate the sensitivity and specificity of the methodology to AD. Global NC was evaluated with a general linear model (GLM) framework by regressing whole-head EEG power envelopes in three frequency bands (theta, alpha and beta) with average fNIRS oxy- and deoxy-hemoglobin concentration changes in the frontal and prefrontal cortices. NC was lower in AD compared to healthy controls (HC) with significant differences in the linkage of theta and alpha bands with oxy- and deoxy-hemoglobin, respectively (p = 0.028 and p = 0.020). Importantly, standalone EEG and fNIRS metrics did not highlight differences between AD and HC. Furthermore, a multivariate data-driven analysis of NC between the three frequency bands and the two hemoglobin species delivered a cross-validated classification performance of AD and HC with an Area Under the Curve, AUC = 0.905 (p = 2.17 × 10−5). The findings demonstrate that EEG-fNIRS may indeed represent a powerful ecological tool for clinical evaluation of NC and early identification of AD.


Introduction
Alzheimer's disease (AD) is a form of dementia that usually begins with slight memory failures and slowly progresses into evident cognitive impairment [1][2][3]. AD is the cause of about 70% of cases of elderly dementia [4]. The underlying mechanisms of the clinical symptoms of AD are poorly understood [5]. However, AD is known to be associated with extracellular amyloid beta deposits [6,7], tau protein abnormalities [8,9] and neuronal loss [10,11]. Moreover, AD has been recently linked to early neurovascular dysfunction, which clearly contributes to disease pathogenesis [12][13][14]. In this perspective, probing functional hyperemia, i.e., neurovascular coupling (NC), could provide relevant diagnostic and prognostic neuromarkers of the disease [15][16][17].
Neurovascular coupling is a physiological mechanism that refers to the causal relationship between local neural activity and subsequent increases in cerebral blood flow (CBF) digit encoding retrieval task. They obtained sufficiently good classification performances (around 80% in a four-class classification) which were superior to those obtained employing standalone modalities. Notably, while they did combine EEG and fNIRS features in a multivariate data-driven analysis, they did not evaluate cross-modality metrics related to NC. In fact, although EEG-fNIRS systems have been employed to extract cross-modality metrics of NC in different clinical populations, such as stroke patients [53], neonates [54], patients with acute brain injury [55] and in the assessment of physiological aging [56], to our knowledge, such multimodal methodology was not applied for the evaluation of NC and its modification in AD.
In this study, a comparison of resting-state EEG-fNIRS NC between AD and healthy controls (HC) was performed. The extent of NC was assessed by correlating whole-head EEG power envelopes in three frequency bands (theta (θ), alpha (α) and beta (β)) with changes in HbO and HbR, assessed through fNIRS in the frontal and prefrontal cortices (six NC metrics). The coupling between electrical and hemodynamic brain activity was evaluated employing a general linear model (GLM) [57][58][59][60]. The GLM is a standard approach to infer brain activity from functional neuroimaging technologies such as fNIRS and fMRI, and it is a generalization of univariate regression. In this context, GLM allowed us to concurrently evaluate the coupling between multiple EEG frequency bands and hemoglobin species. Finally, a multivariate data-driven (i.e., machine learning) analysis of the six global NC metrics was performed to demonstrate the statistical robustness of the findings and to deliver an approach that could support AD diagnosis in clinical settings.

Participants
Thirty-five participants were recruited for the study. The population was composed of 17 AD patients and 18 healthy controls (HC, refer to Table 1 for demographic information). The AD patients had a diagnosis of mild probable Alzheimer's disease, as defined by the Diagnostic and Statistical Manual of Mental Disorders, 5th edition (DSM-5). Patients with moderate cognitive impairment (Mini Mental State Examination, MMSE, < 20/30) [61], vascular dementia, behavioral or psychiatric disorders, brain lesions or with a history of stroke or traumatic brain injury were excluded. Table 1. Demographic information of Alzheimer's disease (AD) and healthy controls (HC) involved in the study.

Experimental Design
Resting-state recordings of 5 min in duration were performed after the administration of the MMSE. The patients sat comfortably on a chair in front of a doctor, and they were asked to stay as still as possible with their eyes closed and to not think anything specific during the recordings.

EEG Instrumentation and Measurements
A high-density, 128-channel, full-head EEG system (Electrical Geodesic Inc., Eugene, OR, USA, EEG System Net 300, Figure 1a) was used to record brain electrical activity. Skin/electrode impedance was measured before recordings and kept below 50 kΩ. Notably, 50 kΩ is considered an impedance within the optimal range for the EEG system employed. Although a sensor-to-scalp impedance below 5 kΩ is generally required for EEG recordings, the HydroCel Geodesics Sensor Net provides excellent signals with impedances up to 50-100 kΩ thanks to the high input impedance amplifiers [62]. EEG data were sampled at 250 Hz.

FNIRS Instrumentation and Measurements
A frequency domain near-infrared spectroscopy (NIRS) system (Imagent, ISS Inc., Champaign, IL, USA) was used to perform the optical measurements. The system was made of 32 laser diodes sources, 16 emitting at 690 nm wavelength and 16 emitting at 830 nm wavelength, and 4 photomultiplier tube (PMT) detectors. The lasers were modulated at 110 MHz, whereas the PMTs were modulated at 110 MHz and 5 KHz for heterodyne detection of the modulated light. A time multiplexing scheme was used on the sources to avoid their crosstalk. The measurement sampling rate (accounting for source time-multiplexing) was set at 10.42 Hz. The light intensity delivered to the scalp was below 4 mW/cm 2 , within the ANSI standard safe limits. NIR light was carried to the scalp using optic fibers with a 400 µm core (2 fibers reached each light injection location delivering light from one laser diode emitting at 830 nm wavelength and one laser diode emitting at 690 nm wavelength) and from the scalp back to the PMTs using fiber bundles with a diameter of 3 mm.
The fibers were located on the frontal and prefrontal areas and were held in place with a home-made optical patch located on top of the high-density EEG cap (Figure 1a). The fiber heads reached the scalp using the space among the electrodes of the EEG cap, that were also employed for positioning of the optical array with reference to the 10/20 system [63] ( Figure 1b). Notably, the need to develop a home-made optical patch was driven by the necessity to locate fNIRS optodes in between EEG electrodes that were arranged in a high-density configuration and that were around 1 cm thick. These mechanical constraints led us to develop a home-made patch that was able to hold the optodes on top of the EEG cap while allowing them to reach the scalp surface for a correct optode-to-scalp coupling using the small space available among EEG electrodes. The optical array allowed us to perform fNIRS measurements from 16 channels at source-detector distances of 35 mm (long separation channels), probing the brain cortex (Figure 1b,c), and from 4 short separation channels at a 15 mm interoptode distance ( Figure 1b). These short separation channels were used to selectively measure hemoglobin concentration changes in the scalp and to correct for their contamination in the long separation channels [64,65].
Slightly longer interoptode distances were employed compared to standard practices for both long and short separation channels (35 mm vs.~30 mm and 15 mm vs. 10 mm) [64]. These interoptode distances were preferred because of the age range and the disease of the population examined. Brain atrophy associated with aging and AD tends to increase the distance of the cortical surface from the scalp. Moreover, only 4 short separation channels were used because of the limited number of optodes available. However, because of the not exceedingly large area covered by the optical array, the average distance between the centroid of the closest short separation channel and the centroid of each long separation channel was 18 mm, which is considered a good compromise [66]. the centroid of the closest short separation channel and the centroid of each long separation channel was 18 mm, which is considered a good compromise [66].

Figure 1. (a)
Electroencephalography and functional near-infrared spectroscopy (EEG-fNIRS) probes located on a participant's head. The fNIRS fiber heads reached the scalp using the spaces between the electrodes of the EEG cap. (b) EEG-fNIRS probes overlaid on a template head. The high-density (128-channel) EEG layout was arranged in agreement with the 10/20 system and the fNIRS probes were positioned with reference to the EEG sensors. (c) Channels' average sensitivity map of the fNIRS measurements. The sensitivity map was computed using the finite element method to solve the diffusion equation [67] and it is displayed up to an attenuation of 60 dB from its maximum value.

EEG Pre-Processing
The flowchart of EEG pre-processing is reported in Figure 2a. EEG data were visually inspected to reject saturated or corrupted epochs. Data were then band-pass filtered between 1 and 80 Hz and notch-filtered at 50 Hz (zero-lag 2nd order Butterworth digital filters). Furthermore, a semiautomatic procedure based on independent component analysis (ICA) was applied to identify and to remove cardiac and ocular artifacts, as well as signals coming from muscles [68,69]. The filtered and artifact-corrected EEG signals were decomposed in three frequency bands of interest (i.e., theta: 3.5-8.2 Hz, alpha: 7.4-13 Hz and beta: 13-30 Hz) and the power temporal envelopes were computed as the absolute values of their Hilbert transform.

FNIRS Pre-Processing
The flowchart of fNIRS pre-processing is reported in Figure 2b. The raw continuouswave components of the acquired signals were converted into optical densities (ODs) according to the equation: where I(t) is the time dependence of the recorded signal intensity and Im is its average value during rest. Furthermore, motion artifacts were corrected by applying a wavelet- Figure 1. (a) Electroencephalography and functional near-infrared spectroscopy (EEG-fNIRS) probes located on a participant's head. The fNIRS fiber heads reached the scalp using the spaces between the electrodes of the EEG cap. (b) EEG-fNIRS probes overlaid on a template head. The high-density (128-channel) EEG layout was arranged in agreement with the 10/20 system and the fNIRS probes were positioned with reference to the EEG sensors. (c) Channels' average sensitivity map of the fNIRS measurements. The sensitivity map was computed using the finite element method to solve the diffusion equation [67] and it is displayed up to an attenuation of 60 dB from its maximum value.

EEG Pre-Processing
The flowchart of EEG pre-processing is reported in Figure 2a. EEG data were visually inspected to reject saturated or corrupted epochs. Data were then band-pass filtered between 1 and 80 Hz and notch-filtered at 50 Hz (zero-lag 2nd order Butterworth digital filters). Furthermore, a semiautomatic procedure based on independent component analysis (ICA) was applied to identify and to remove cardiac and ocular artifacts, as well as signals coming from muscles [68,69]. The filtered and artifact-corrected EEG signals were decomposed in three frequency bands of interest (i.e., theta: 3.5-8.2 Hz, alpha: 7.4-13 Hz and beta: 13-30 Hz) and the power temporal envelopes were computed as the absolute values of their Hilbert transform. The end results of the EEG pre-processing were artifact-corrected EEG power envelopes for the 128 channels in the three frequency bands of interest (theta, alpha and beta). (b) Pre-processing of the fNIRS involving Optical Densities (OD) computation, wavelet-based motion correction and oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) computation. The end results of the fNIRS pre-processing were artifactcorrected HbO and HbR changes in both the long and the short separation channels.

Neurovascular Coupling Estimation
The algorithmic approach for NC estimation is reported in Figure 3. For each subject and frequency band, EEG power timecourses were normalized (z-scored) and fed to a temporal principal component analysis (PCA). This analysis provided temporally independent timecourses (principal components, PCs) that were common to the different EEG channels within each band. The first PC, i.e., the component that explained the largest Figure 2. (a) Pre-processing of the EEG involving band-pass filtering, independent component analysis (ICA) motion correction and power envelope computation. The end results of the EEG pre-processing were artifact-corrected EEG power envelopes for the 128 channels in the three frequency bands of interest (theta, alpha and beta). (b) Pre-processing of the fNIRS involving Optical Densities (OD) computation, wavelet-based motion correction and oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) computation. The end results of the fNIRS pre-processing were artifact-corrected HbO and HbR changes in both the long and the short separation channels.

FNIRS Pre-Processing
The flowchart of fNIRS pre-processing is reported in Figure 2b. The raw continuouswave components of the acquired signals were converted into optical densities (ODs) according to the equation: where I(t) is the time dependence of the recorded signal intensity and I m is its average value during rest. Furthermore, motion artifacts were corrected by applying a wavelet-based procedure [70] and the ODs were band-pass filtered with a zero-lag, 4th order Butterworth digital filter with cut-off frequencies of 0.01 Hz and 0.4 Hz. These motion correction and filtering approaches allowed us to highlight the hemodynamic components of the acquired optical signal. Notably, the recordings were performed standing still and at rest and they were poorly contaminated by motion-related signals. Thus, the use of a motion correction algorithm that does not distort the data in the presence of a low level of motion noise was preferred (kurtosis-based wavelet filtering) [70]. After visual inspection of the filtered and motion-corrected ODs, changes in the concentration of HbO and HbR were derived for the rest period from each channel based on the modified Lambert-Beer law [71]: where ρ is the interoptode distance (either 35 mm or 15 mm), ε is the extinction coefficient for each chromophore and wavelength and DPF is the differential pathlength factor for each wavelength. The extinction coefficients were extracted from Zijlstra et al. [72], whereas the DPFs [73], accounting not only for wavelength used but also for subject's age, were extracted from Scholkmann and Wolf [74].

Neurovascular Coupling Estimation
The algorithmic approach for NC estimation is reported in Figure 3. For each subject and frequency band, EEG power timecourses were normalized (z-scored) and fed to a temporal principal component analysis (PCA). This analysis provided temporally independent timecourses (principal components, PCs) that were common to the different EEG channels within each band. The first PC, i.e., the component that explained the largest amount of signal variance, was retained for each frequency band. This analysis effectively provided a global EEG signal for each frequency band which was not obtained through simple arithmetic averaging among channels but through weighted averaging that accounted for temporal correlation among the different channels. As the EEG signal is highly correlated among channels due to volume conduction effects, this approach is intrinsically robust to non-homogenous noise across channels.
The first PCs of the three EEG frequency bands were convolved with a canonical HRF [24] to obtain a model of the expected hemoglobin concentration changes following resting-state electrical brain activity. These three global signals were resampled (at 1 Hz) and normalized (z-score) and then they were considered as independent variables of the GLM. The normalized (z-score) and resampled (at 1 Hz) HbO and HbR changes in the long separation fNIRS channels were employed as dependent variables. In order to regress out the scalp contamination in the long separation channels, the normalized hemoglobin changes in the four short separation fNIRS channels were also considered as independent variables. The NC estimate was encoded in the β-weight, that, given the normalization of the signals prior to the GLM, delivered a standardized regression coefficient of each EEG power envelope association with the HbO and HbR changes for each of the long separation fNIRS channels. Notice that, whereas the GLM has already been used in EEG-fMRI studies [75,76], it is novel in the EEG-fNIRS investigation of brain activity and NC. tion fNIRS channels. Notice that, whereas the GLM has already been used in EEG-fMRI studies [75,76], it is novel in the EEG-fNIRS investigation of brain activity and NC.
In order to prove the utility of the EEG-fNIRS multimodal approach and GLM for studying the NC in AD, unimodal electrical and hemodynamic brain activities and their changes with AD were investigated as well. In particular, the average EEG and fNIRS powers were evaluated by measuring the mean value of the EEG power in the theta, alpha and beta bands and the standard deviation in the HbO and HbR fNIRS timecourses.  In order to prove the utility of the EEG-fNIRS multimodal approach and GLM for studying the NC in AD, unimodal electrical and hemodynamic brain activities and their changes with AD were investigated as well. In particular, the average EEG and fNIRS powers were evaluated by measuring the mean value of the EEG power in the theta, alpha and beta bands and the standard deviation in the HbO and HbR fNIRS timecourses.

Statistical Inference and Multivariate Data-Driven Analysis
The global NC of the three EEG frequency bands with HbO and HbR (6 metrics of NC for each subject) were evaluated by computing the average of the GLM β-weights among the fNIRS long separation channels.
The global electrical and hemodynamic brain activities were estimated by computing the channels' average EEG power in the theta, alpha and beta bands (3 metrics of electrical brain activity for each subject) and the channels' average fNIRS power (i.e., signal standard deviation) for HbO and HbR timecourses (2 metrics of hemodynamic brain activity for each subject).
In order to test the statistical significance of NC in AD and HC and test between-groups differences, t-tests were performed [77].
Once univariate statistical significance was explored, in order to further corroborate the results and to explore the clinical significance of the findings, a multivariate data-driven analysis was implemented [78]. A multivariate regression analysis was performed in the attempt to classify each subject as AD or HC based on EEG, fNIRS and NC metrics. Since the classification was known, such an approach represented a classical problem of supervised learning in a machine learning framework [79]. Although several machine learning approaches could be suitable for such a purpose, given the limited number of independent features available and the exploratory nature of the implemented approach, the procedural complexity was limited to the use of a simple linear regressor followed by an output threshold selection. Moreover, the use of a linear regression avoided investigation of the model hyperparameter space, which is compulsory if regularization approaches or non-linear regressors are used (e.g., choosing the regularization parameter or the nonlinearity extent). The closed solution and the small number of effective degrees of freedom of the linear regressor (equal to the number of the independent features) greatly reduced the chance of overfitting [80].
In order to build and test the multivariate linear regression, a value of "1" and "0" was attributed to AD and HC, respectively, and the model was trained on all subjects except one and tested on the remaining subject in an iterative manner (i.e., using the leave-one-out cross-validation (CV) scheme) [81,82]. The classification outcome was assessed through a receiver operating characteristic (ROC) analysis of the out-of-training-sample outputs of the iteratively trained regressor which, given the labeling value associated with AD and HC, tended to provide larger values for the AD group. The framework was implemented multiple times to estimate the capability of the procedure to separate AD from HC using the 6 global EEG-fNIRS NCs, the 3 global EEG powers in the theta, alpha and beta bands and the 2 global fNIRS powers for the HbO and HbR timecourses.
In order to investigate the effect of the EEG channel distribution on the global characteristics of the NC metrics, the exact same multivariate data-driven analysis using NC metrics was performed again, but selectively evaluating the EEG signals acquired with half of the available electrodes, all located on the frontal and prefrontal brain areas, while excluding the occipital ones.
Finally, in order to evaluate the performance of the proposed approach as a function of practical variables of interest, such as EEG-fNIRS recording time and number of EEG channels employed, the multivariate data-driven analysis was conducted multiple times by changing either the temporal duration of the signals used (from 15 s to 300 s with a 15 s step) or the number of channels employed (from 4 to 128 with a 4-channel step). The different time windows were selected at the center of the recordings for each subject, whereas the subsample of channels was randomly selected from the 128 channels covering the whole head. When the 300 s window or the 128 channels were evaluated, the analyses were equivalent to the original one where all the signals recorded were employed for NC evaluation and AD vs. HC classification. Table 1 reports demographic information for AD and HC. Notably, AD and HC were matched for gender, age and education. A lower MMSE score was identified in AD (t = 5.51, df = 33, p = 4.05 × 10 −6 ).

Results
The first PCs of the EEG power envelopes explained an average 64% (SD = 13%), 68% (SD = 14%) and 69% (SD = 12%) of the variance in the theta, alpha and beta bands, respectively. Figure 4 shows the distributions of the GLM standardized β-weights for the long separation fNIRS channels in the different subjects (reported with different colors) for the three EEG frequency bands and HbO and HbR in AD and HC.
The first PCs of the EEG power envelopes explained an average 64% (SD = 13%), 68% (SD = 14%) and 69% (SD = 12%) of the variance in the theta, alpha and beta bands, respectively. Figure 4 shows the distributions of the GLM standardized β-weights for the long separation fNIRS channels in the different subjects (reported with different colors) for the three EEG frequency bands and HbO and HbR in AD and HC.   Figure 5a reports the NC of theta, alpha and beta bands with HbO for AD and HC, whereas Figure 5b reports the NC of theta, alpha and beta bands with HbR for AD and HC.
A statistically significant NC was obtained for HC in the coupling between theta band and HbO (t = 2.12, df = 17, p = 0.049) and in the coupling between the alpha band and HbR (t = 2.24, df = 17, p = 0.039). This link was lost in AD (p > 0.05). In fact, a statistically significant decrease in global NC in AD with respect to HC was obtained for the coupling between the theta band and HbO (t = −2.30, df = 33, p = 0.028) and for the alpha band and HbR (t = −2.44, df = 33, p = 0.020). Notably, no NC was statistically above zero and no significant increase in global NC was obtained for AD with respect to HC.   Figure 5a reports the NC of theta, alpha and beta bands with HbO for AD and HC, whereas Figure 5b reports the NC of theta, alpha and beta bands with HbR for AD and HC. For comparison, Figure 6 reports the boxplots showing global metrics of electrical and hemodynamic brain activity evaluated through standalone EEG and fNIRS. Figure 6a reports the global EEG power in theta, alpha and beta bands for AD and HC. Figure 6b A statistically significant NC was obtained for HC in the coupling between theta band and HbO (t = 2.12, df = 17, p = 0.049) and in the coupling between the alpha band and HbR (t = 2.24, df = 17, p = 0.039). This link was lost in AD (p > 0.05). In fact, a statistically significant decrease in global NC in AD with respect to HC was obtained for the coupling between the theta band and HbO (t = −2.30, df = 33, p = 0.028) and for the alpha band and HbR (t = −2.44, df = 33, p = 0.020). Notably, no NC was statistically above zero and no significant increase in global NC was obtained for AD with respect to HC.
For comparison, Figure 6 reports the boxplots showing global metrics of electrical and hemodynamic brain activity evaluated through standalone EEG and fNIRS. Figure 6a reports the global EEG power in theta, alpha and beta bands for AD and HC. Figure 6b reports the global HbO and HbR power for AD and HC. No statistically significant differences were found between AD and HC for the metrics evaluated through unimodal recordings.  Figure 7 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the six NC metrics. Importantly, the results were obtained using a leave-one-out CV, depicting the generalization performance of the approach. Figure 7a reports the outcome of the ROC analysis where an area under the curve (AUC) of 0.905 (z = 4.090, p = 2.17 × 10 −5 ) was obtained. Table 2 reports the confusion matrix that was computed by selecting the threshold of the outcome of the multivariate linear regression that provided a sensitivity of 88.2% and a specificity of 94.5%. Figure 7b shows the CV loadings of the multivariate regression analysis, showing how the different features were concurrently weighted to deliver the classification. Since in the classification procedure AD had a value of "1", whereas HC had a value of "0", a positive load depicts a larger NC for AD whereas a negative load represents a smaller NC for AD. Notice that, whereas the upper panel of Figure 7b reports the regression loadings, the lower panel of Figure 7b reports the loadings normalized by their variability (i.e., their standard deviation) estimated using a bootstrap approach (10,000 iterations) and effectively providing the statistical significance (z-score value) of each multivariate load.  Figure 7 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the six NC metrics. Importantly, the results were obtained using a leaveone-out CV, depicting the generalization performance of the approach. Figure 7a reports the outcome of the ROC analysis where an area under the curve (AUC) of 0.905 (z = 4.090, p = 2.17 × 10 −5 ) was obtained. Table 2 reports the confusion matrix that was computed by selecting the threshold of the outcome of the multivariate linear regression that provided a sensitivity of 88.2% and a specificity of 94.5%. Figure 7b shows the CV loadings of the multivariate regression analysis, showing how the different features were concurrently weighted to deliver the classification. Since in the classification procedure AD had a value of "1", whereas HC had a value of "0", a positive load depicts a larger NC for AD whereas a negative load represents a smaller NC for AD. Notice that, whereas the upper panel of Figure 7b reports the regression loadings, the lower panel of Figure 7b reports the loadings normalized by their variability (i.e., their standard deviation) estimated using a bootstrap approach (10,000 iterations) and effectively providing the statistical significance (z-score value) of each multivariate load.  For comparison with the multivariate data-driven analysis of NC, Figure 8 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the standalone EEG and fNIRS metrics. Figure 8a reports the CV outcome of the ROC analysis employing the three EEG global powers, whereas Figure 8b reports the results of the ROC analysis employing the two fNIRS global hemoglobin powers. Using unimodal analysis, a poor AD vs. HC classification was obtained with AUC = 0.521 (z = 0.212, p = 0.416) and AUC = 0.555 (z = 0.555, p = 0.289) for the EEG and the fNIRS features, respectively. Importantly, the multivariate classification performance based on the multimodal NC evaluation was statistically superior to both the unimodal EEG and fNIRS analysis (AUC = 0.905 vs. 0.521, z = 3.403, p = 6.6 × 10 −4 ; AUC = 0.905 vs. 0.555, z = 3.115, p = 1.8 × 10 −3 ).  For comparison with the multivariate data-driven analysis of NC, Figure 8 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the standalone EEG and fNIRS metrics. Figure 8a reports the CV outcome of the ROC analysis employing the three EEG global powers, whereas Figure 8b reports the results of the ROC analysis employing the two fNIRS global hemoglobin powers. Using unimodal analysis, a poor AD vs. HC classification was obtained with AUC = 0.521 (z = 0.212, p = 0.416) and AUC = 0.555 (z = 0.555, p = 0.289) for the EEG and the fNIRS features, respectively. Importantly, the multivariate classification performance based on the multimodal NC evaluation was statistically superior to both the unimodal EEG and fNIRS analysis (AUC = 0. Notably, when the whole analysis was re-run using half of the total electrodes, eliminating the occipital ones, the difference between AD and HC in the coupling of the theta band and HbO was still statistically significant (t = −2.170, df = 33, p = 0.037), whereas the difference in the coupling between the alpha band and HbR vanished (t = −0.530, df = 33, p = 0.600), with the multivariate data-driven analysis delivering an AUC of 0.830 (z = 3.329 p = 4.4 × 10 −4 ). Finally, Figure 9 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the six NC metrics as a function of either the recording time or the number of randomly selected EEG channels. Figure 9a shows the loadings of the multivariate regression analysis and the associated z-scores as a function of recording time available. Figure 9b instead shows the AUC of the classification outcome, always as a function of the recording time. Indeed, the loadings and the associated z-scores converged to negative values, whereas the AUC increased as a function of recording time. Figure 9c shows the loadings and the associated z-scores as a function of the number of randomly selected EEG channels. Figure 9d shows the AUC of the associated classification outcome. In this case, the loadings quickly assumed negative values, whereas the AUC was quite stable. Notably, when the whole analysis was re-run using half of the total electrodes, eliminating the occipital ones, the difference between AD and HC in the coupling of the theta band and HbO was still statistically significant (t = −2.170, df = 33, p = 0.037), whereas the difference in the coupling between the alpha band and HbR vanished (t = −0.530, df = 33, p = 0.600), with the multivariate data-driven analysis delivering an AUC of 0.830 (z = 3.329 p = 4.4 × 10 −4 ). Finally, Figure 9 reports the outcome of the multivariate data-driven analysis trying to predict AD and HC from the six NC metrics as a function of either the recording time or the number of randomly selected EEG channels. Figure 9a shows the loadings of the multivariate regression analysis and the associated z-scores as a function of recording time available. Figure 9b instead shows the AUC of the classification outcome, always as a function of the recording time. Indeed, the loadings and the associated z-scores converged to negative values, whereas the AUC increased as a function of recording time. Figure 9c shows the loadings and the associated z-scores as a function of the number of randomly selected EEG channels. Figure 9d shows the AUC of the associated classification outcome. In this case, the loadings quickly assumed negative values, whereas the AUC was quite stable.

Identification of EEG-fNIRS De-Coupling in AD
AD is characterized by a progressive cognitive decline over the course of the disease that slowly disrupts daily life [83]. The exact mechanisms inducing AD symptoms remain largely unknown [84]. Nonetheless, novel neuroimaging modalities and advanced data analysis may enhance the ability to non-invasively assess brain status and subsequently identify its functional deterioration in AD. Such approaches must be suitable for ambulatory settings in order to impact standard clinical practice. To this end, multimodal evaluation of brain activity and NC through synchronous EEG-fNIRS recordings represents a promising tool [46]. EEG and fNIRS are relatively cheap, portable and easy to use technologies. These characteristics make them highly suited for both screening-like and longitudinal ambulatory evaluation of brain diseases. Synchronous EEG and fNIRS recordings may provide relevant information regarding NC, whose deterioration is known to be associated with AD [12][13][14].
In this study, multimodal EEG-fNIRS recordings were performed at rest and they were used to infer NC and evaluate its impairment in AD. The multivariate correlations of whole-head EEG power envelopes in the theta, alpha and beta bands with fNIRS-derived changes in HbO and HbR in frontal and prefrontal cortices were estimated by means of a GLM approach and representative metrics of NC. A significant decrease in correlation, i.e., a global neurovascular un-coupling, was found in AD patients between the theta band power envelope and HbO changes (t = −2.30, df = 33, p = 0.028) and between the alpha band power envelope and HbR changes (t = −2.44, df = 33, p = 0.020, Figure 5). Hence, for simplicity, we preferred to only report the results about theta and beta bands in the paper as representative of low-and high-frequency EEG band associations with fNIRS hemoglobin oscillations.

Identification of EEG-fNIRS De-Coupling in AD
AD is characterized by a progressive cognitive decline over the course of the disease that slowly disrupts daily life [83]. The exact mechanisms inducing AD symptoms remain largely unknown [84]. Nonetheless, novel neuroimaging modalities and advanced data analysis may enhance the ability to non-invasively assess brain status and subsequently identify its functional deterioration in AD. Such approaches must be suitable for ambulatory settings in order to impact standard clinical practice. To this end, multimodal evaluation of brain activity and NC through synchronous EEG-fNIRS recordings represents a promising tool [46]. EEG and fNIRS are relatively cheap, portable and easy to use technologies. These characteristics make them highly suited for both screening-like and longitudinal ambulatory evaluation of brain diseases. Synchronous EEG and fNIRS recordings may provide relevant information regarding NC, whose deterioration is known to be associated with AD [12][13][14].
In this study, multimodal EEG-fNIRS recordings were performed at rest and they were used to infer NC and evaluate its impairment in AD. The multivariate correlations of whole-head EEG power envelopes in the theta, alpha and beta bands with fNIRS-derived changes in HbO and HbR in frontal and prefrontal cortices were estimated by means of a GLM approach and representative metrics of NC. A significant decrease in correlation, i.e., a global neurovascular un-coupling, was found in AD patients between the theta band power envelope and HbO changes (t = −2.30, df = 33, p = 0.028) and between the alpha band power envelope and HbR changes (t = −2.44, df = 33, p = 0.020, Figure 5). Hence, for simplicity, we preferred to only report the results about theta and beta bands in the paper as representative of low-and high-frequency EEG band associations with fNIRS hemoglobin oscillations.
The selective un-coupling between specific EEG frequency bands and fNIRS hemoglobin species requires further discussion. The decrease in EEG-fNIRS correlations in AD was found for EEG frequency bands and hemoglobin species where the average NC was significantly above chance in HC (i.e., average β-weight larger than zero, Figure 5). Notably, these results are in line with previous multimodal works combining EEG with BOLD fMRI in HC [45,85].
In HC, a strong negative correlation of EEG alpha power with parietal and frontal cortical activity of BOLD was found [86], whereas a positive correlation was identified at rest between BOLD and theta power of local field potentials in parahippocampal areas [87]. Since the BOLD signal has a coarse negative association with HbR and positive association with HbO [23], both the fMRI findings and the results of this study suggest that the selective de-coupling found in AD between theta and HbO and alpha and HbR could simply reflect a general neurovascular un-coupling that is more pronounced for frequency bands and hemoglobin species where the original physiological association is stronger.
Moreover, none of the six EEG-fNIRS NC metrics were significantly higher in AD with respect to HC. Indeed, these results supported the hypothesis of NC decline in AD. To further highlight the importance of investigating NC in AD with multimodal EEG-fNIRS recordings, unimodal evaluation of either global electrical or hemodynamic brain activity, expressed as average EEG and fNIRS signal powers, did not highlight statistically significant differences between AD and HC. Although not statistically significant, it is interesting to notice that a qualitatively lower EEG power in AD with respect to HC was found for all frequency bands considered coupled with a qualitatively higher fNIRS power for both HbO and HbR. This result may be indicative of a dysregulation of NC which is selectively identified by the multimodal analysis.

Multivariate Exploitation of EEG-fNIRS De-Coupling in AD
To further investigate the NC dysregulation in AD and to exploit the potentialities of the approach presented, a multivariate data-driven analysis was also performed. This analysis concurrently evaluated how the NC metrics (i.e., the standardized regression weights between the three frequency bands and the two hemoglobin species, six features) were modified in AD and it investigated to what extent these NC metrics could be used to identify the pathology. The procedure classified AD and HC relying on the six NC metrics and a multiple linear regression approach. The CV classification outcome was assessed by means of an ROC analysis that delivered AUC = 0.905 (z = 4.090, p = 2.17 × 10 −5 , Figure 7a). An interesting result of the multivariate data-driven analysis was that the expected value of the loadings of the linear regression, associated to each of the six NC features was always negative (Figure 7b). This finding further supported the neurovascular un-coupling in AD already found in the univariate analysis. Moreover, the highest statistical relevance of the magnitude of the multivariate loadings was indeed found for the two NCs that were significantly lower in AD based on the univariate analysis (i.e., the link between the theta band and HbO and the link between the alpha band and HbR).
Always in accordance with the univariate analysis, the multivariate data-driven approach used on standalone EEG and fNIRS powers did not provide a classification accuracy of AD and HC above chance (Figure 8). These results further demonstrated the higher sensitivity and specificity to AD of NC metrics evaluated through multimodal approaches with respect to brain activity metrics estimated based on unimodal procedures (AUC = 0.905 vs. 0.521, z = 3.403, p = 6.6 × 10 −4 ; AUC = 0.905 vs. 0.555, z = 3.115, p = 1.8 × 10 −3 ). At a specific threshold of the outcome of the multivariate data-driven analysis of NC, a sensitivity of 88.2% and a specificity of 94.5% were obtained, providing additional evidence of the high potentialities of the method.
Of note, whereas the univariate evaluation of the couplings between multiple EEG frequency bands and fNIRS hemoglobin forms might have produced false positives, the multivariate analysis indeed relied on considering all the information at once, greatly decreasing the chance of false positives related to multiple comparison. Indeed, the highly significant multivariate results support the positive findings of the study.
As a further investigation of the findings, the multivariate classification of AD and HC based on EEG-fNIRS NC metrics was evaluated as a function of either the recording time available or the number of EEG channels employed (Figure 9). These are important aspects of the replicability of the results and the implementation of the approach in a clinical environment. With respect to recording time length, both the multivariate loadings and associated z-scores tended to reach stable and negative values after 200 s. The AUC was close to the best performance only when the recording time was close to 300 s, while showing a rather linear decrease in value with shorter recording times. These results corroborated the robustness of the original findings but also highlighted the need for recording durations of several minutes to employ the proposed approach with sufficient sensitivity and specificity. On the contrary, both the multivariate loadings and associated z-scores were stable with the number of randomly selected EEG channels and, importantly, the AUC was very close to the best classification performance even with particularly small numbers of channels used. These findings highlighted the applicability of the approach in a clinical environment where a small number of electrodes is generally available.

Study Limitations
Several limitations of the study must be commented upon and might be addressed in future studies.
Firstly, whereas the high-density EEG provided whole-head coverage, the fNIRS layout allowed us to limit the investigation to the frontal and prefrontal cortices. The layout choice was driven by practical limitations of the available fNIRS instrumentation that provided a limited number of optodes. The choice to cover the frontal and prefrontal cortices was led by the higher sensitivity of the fNIRS optodes in areas not covered by hair and by the involvement of these areas in AD-related cognitive impairment, although to a lower extent with respect to other areas, such as the medial temporal lobe [88,89]. In order to control for the different head coverage levels of EEG and fNIRS, the analysis was run again, but by selectively choosing EEG signals acquired with half of the electrodes all located on the frontal and prefrontal brain areas. With this analysis, the difference between AD and HC in the coupling of the theta band and HbO remained unaltered (t = −2.170, df = 33, p = 0.037), whereas the difference in the coupling between the alpha band and HbR still pointed toward a reduction of the linkage strength in AD but vanished statistically (t = −0.530, df = 33, p = 0.600). Because of the reduced effect in the coupling between the alpha band and HbR, the multivariate data-driven analysis delivered an AUC of 0.830, which was below the original one of 0.905. These findings suggest that the analysis of the whole-head EEG signal identified an increase in global alpha activity (known to be associated with brain electrical activity de-activation) that induced a decrease in HbR (known to be associated with brain hemodynamic de-activation) in frontal and prefrontal areas and that the alteration of this effect in AD was useful in identifying the patients with the disease. This information was lost when limiting the useful electrodes to those located in the frontal region because of the diminished sensitivity of those electrodes to the alpha band activity. This result, coupled with the findings regarding the classification performance as a function of the number of randomly selected EEG channels, suggested that, although a limited number of electrodes can be used to exploit the proposed approach, such electrodes should be evenly distributed across the scalp. However, further studies should be performed with a whole-head fNIRS system to assess a global NC estimated using whole-brain hemodynamic information in conjunction with electrical activity. Nonetheless, it might be possible that a global NC assessed through whole-head fNIRS may deliver a lower sensitivity for AD. In fact, since the pathology is associated with a cognitive decline often related to frontal cortex activity dysregulation, AD may be characterized by a selective disruption of fNIRS modulation within the frontal lobes in response to global EEG oscillations [90]. Particularly, it could be worth extending the fNIRS investigation to other brain regions highly impaired by AD, such as the medial temporal lobe [89].
Secondly, always with reference to the electrode and optode distribution on the scalp, a less relevant issue, mainly driven by mechanical limitations, was associated with the fact that the EEG electrodes where not placed on the same scalp locations as the optodes nor were exactly in between an optical source and a detector. However, we should further stress that EEG, due to volume conduction effects, is intrinsically a low spatial resolution technology. Indeed, a displacement of EEG sensors of a few centimeters does not drastically modify the signal recorded. In this perspective, a small displacement between fNIRS optodes and EEG electrodes does not affect the assessment of a global NC using EEG-fNIRS, also considering that the global signal was extracted using a high-density EEG system which has an average distance between sensors well below the minimum wavelength (corresponding to the maximum spatial frequency component) of the measured voltages across the head [91]. This statement was indeed heavily supported by the results regarding the stability of the proposed approach as a function of the number of randomly selected EEG channels used.
Thirdly, although the optical probes were located with reference to the EEG 10/20 system, an accurate co-registration of the optical sensors with a specific subject anatomy (e.g., structural MRI of the head) was missing. In fact, also for this reason, only global estimates of NC were evaluated, effectively not accounting for its spatial characteristics. Further studies should indeed be focused on assessing the spatial features of NC modifications in AD. However, although the spatial characteristics of NC may provide additional information, a global metric of NC can be more directly employed in a clinical environment. In fact, providing an accurate spatial localization of electrical and hemodynamic brain activity with EEG and fNIRS requires several methodological steps that make such an approach impractical for employment outside research environments (e.g., co-registration with an anatomical image of the subject, implementation of forward and inverse approaches [92,93]).
Fourthly, although the multivariate data-driven analysis delivered a high classification performance, further studies should be conducted to increase the overall sample size of the population. In fact, these approaches rely on analyses that might drastically increase their performances when large sample sizes are available [80]. When a larger data sample is available, the analysis of NC might be extended beyond the six global metrics that were used in this work, for example, considering more than one principal component for each EEG frequency band, and with the possibilities to explore non-linearities in the data.
Fifthly, in vivo neuroimaging metrics, such as the NC measured with multimodal EEG-fNIRS recordings, are certainly not suitable to understand the underlying microscopic mechanisms associated with the macroscopic manifestations. The drop in correlation between EEG and fNIRS signals identified in AD patients was plausibly associated with a disruption of signal pathways in astrocytes, vascular smooth muscle cells and endothelia, causing a decreased or slower cerebrovascular reactivity in NC, but further study should be performed to fully explain the findings [15].
Finally, this study was limited to the comparison of NC between AD and HC. In order to employ this approach for diagnostic purposes, the capabilities of this method to identify AD and to differentiate it from other similar pathologies should be assessed, effectively providing a procedure able to perform a differential diagnosis [94]. This is particularly relevant since other forms of dementia might indeed exhibit similar neurovascular uncoupling to that found in AD. In fact, since AD shares similar symptomologies with other neurological pathologies, early diagnosis and monitoring of AD represents a major issue for clinicians.

Study Significance
The findings clearly demonstrated a neurovascular de-coupling in AD and highlighted how this de-coupling might be exploited to identify AD at an early stage of the disease. Importantly, the approach implemented can be easily translated to clinical, ambulatory settings, thanks to the portability and usability of both EEG and fNIRS. Thanks to the extraction of global metrics of brain activity to assess NC, this procedure does not require EEG and fNIRS sensor co-registration to an anatomical image, which would imply the usage of expensive and time-consuming imaging techniques (e.g., MRI or CT). The unbiased cross-validated framework built to assess the performance of the multivariate approach in identifying AD through global EEG-fNIRS NC clearly suggests the reliability of the proposed method, despite the small number of subjects included in the study.

Future Perspectives
The neuroimaging approaches presented in this work, unifying novel ecological multimodal resting-state recordings of brain activity and multivariate data-driven analysis of the signals acquired, could support early AD diagnosis, potentially improving treatment efficacy and decreasing societal economic effort. In fact, the actual measurement does not require specifically trained personnel, and the algorithm, once implemented, can be easily used by clinicians.

Conclusions
In this study, NC was evaluated in AD and HC by means of multimodal EEG-fNIRS recordings performed at rest. Six metrics of global NC of theta, alpha and beta band power envelopes with HbO and HbR changes were extracted by using whole-head high-density EEG recordings, frontal and prefrontal cortex fNIRS measures and a GLM approach. Univariate analysis of the six NC features showed a statistically significant decrease in the coupling between theta and HbO and between alpha and HbR. Importantly, no significant effects were found when analyzing the average power of the signals of standalone modalities. A linear, multivariate and data-driven approach to the classification of AD and HC from the six NC metrics delivered a high classification performance with AUC = 0.905. Moreover, the procedure provided negative loadings for all the features, suggesting a general neurovascular un-coupling in AD. Although such an approach should be tested by taking into account other similar pathologies that can affect NC in order to provide a differential diagnosis, the findings indeed demonstrated the capabilities of resting-state EEG-fNIRS to assess NC in ecological ambulatory settings and corroborated the hypothesis of NC disruption in AD, paving the way to the employment of such neuroimaging approaches to support early AD diagnosis.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy issues.

Conflicts of Interest:
The authors declare no conflict of interest.