Comparison of group-level, source localized activity for simultaneous functional near-infrared spectroscopy-magnetoencephalography and simultaneous fNIRS-fMRI during parametric median nerve stimulation

Abstract. Functional near-infrared spectroscopy (fNIRS) is a noninvasive neuroimaging technique, which uses light to measure changes in cerebral blood oxygenation through sensors placed on the surface of the scalp. We recorded concurrent fNIRS with magnetoencephalography (MEG) and functional magnetic resonance imaging (fMRI) in order to investigate the group-level correspondence of these measures with source-localized fNIRS estimates. Healthy participants took part in both a concurrent fNIRS-MEG and fNIRS-fMRI neuroimaging session during two somatosensory stimulation tasks, a blocked design median nerve localizer and parametric pulsed-pair median nerve stimulation using interpulse intervals from 100 to 500 ms. We found the spatial correlation for estimated activation patterns from the somatosensory task was R=0.54, 0.57, and −0.48 and the amplitude correlation was R=0.80, 0.52, and −0.70 for fMRI-MEG, fMRI-fNIRS oxy-hemoglobin, and fMRI-fNIRS deoxy-hemoglobin signals, respectively. Taken together, these results show good correspondence among the fMRI, fNIRS, and MEG with the great majority of the difference across modalities being driven by lower sensitivity for deeper brain sources in MEG and fNIRS. These results provide an important validation of source-localized fNIRS in the context of concurrent multimodal imaging for future studies of the relationship between physiological effects in the human brain.

Comparison of group-level, source localized activity for simultaneous functional near-infrared spectroscopy-magnetoencephalography and simultaneous fNIRS-fMRI during parametric median nerve stimulation Theodore Huppert, a,b, * Jeff Barker, b Benjamin Schmidt, b Shawn Walls, c and Avniel Ghuman c,d

Introduction
Functional near-infrared spectroscopy (fNIRS) is a noninvasive brain imaging technique that uses visible red to near-infrared light (650 to 900 nm) to measure the spectroscopic optical properties of tissue. In particular, for fNIRS brain imaging, light is used to monitor changes in blood oxygenation in the brain through sensors placed on the surface of the scalp. Hemoglobin changes are estimated based on the differential optical absorption profiles of oxygenated and deoxygenated blood (oxy-and deoxy-hemoglobin, HbO 2 , and Hb, respectively). Since the introduction of this method by Jobsis, 1 diffuse optical techniques have been applied to a large variety of different brain-related topics in both clinical and research applications. Several recent reviews have detailed the use of fNIRS brain imaging in the fields of child 2,3 and infant 4,5 research, clinical monitoring, 6,7 and adult brain imaging studies. 8,9 Functional NIRS can play an important role in multimodal imaging because its portability readily allows for simultaneous acquisition with larger, immobile neuroimaging methods such as functional magnetic resonance imaging (fMRI) or magnetoencephalography (MEG). In general, multimodal studies can pair modalities based on the compliment of their particular strengths and weaknesses. MEG acquires high temporal resolution images of the brain's electrophysiological response, but is not sensitive to hemodynamic changes and is restrictive with regards to subject movement. In comparison, fNIRS allows a more portable means to measure cerebral hemodynamic changes and provides additional information about both oxygenated and deoxygenated forms of hemoglobin, but is limited its temporal resolution compared to neural measures. Multimodal brain imaging methods using fNIRS in combination with fMRI, MEG, or electroencephalography (EEG) can lead to a better understanding the physiological basis of human brain activity through investigation of the neural, metabolic, and vascular relationships in the brain. Functional NIRS uses fiber optic cables, which allow remote imaging of a participant inside either fMRI or MEG scanners. Thus, fNIRS can provide a common measurement between sessions to bridge multimodal information. Not only can these methods provide insight into both neural and vascular functionality of the brain, but also can provide independent *Address all correspondence to: Theodore Huppert, E-mail: huppertt@upmc .edu measurements of the signal allowing for technique cross validation in future studies.
While to date, there have been numerous comparisons of concurrent fMRI and fNIRS (reviewed by Steinbrink et al. 10 ), most of these studies have focused on the temporal correlation of these modalities. In general, a consensus of strong linear correlations between the fNIRS and particular deoxy-hemoglobin component and the fMRI-BOLD signal has been observed. However, there has been considerably less work to examine the spatial agreement of these methods beyond qualitative comparisons. In particular, the reconstruction of images of brain activity in fNIRS is an ill-poised and under-determined problem that requires advanced mathematical techniques to simulate the sensitivity (forward) measurement models of the diffusion of light in the head and preform inversion of the model to construct images. In recent years, our group has developed a number of image reconstruction methods based on cortical surface topology for improving image reconstruction of fNIRS data. [11][12][13][14][15] The purpose of this study was to evaluate the performance of these recent image reconstruction methods applied to low-density fNIRS measurements and compare these results to fMRI and MEG source localization. Healthy volunteer participants completed two neuroimaging sessions involving concurrent fNIRS/MEG and concurrent fNIRS/fMRI during a somatosensory stimulation task. Two experimental conditions were examined. First, a 10-s median-nerve localizer task was used to examine the spatial agreement among these three modalities. The spatial location and overlap of the estimate of brain activity from all three modalities was examined in the localizer task. Second, a parametric pulsed-pair study was conducted to look at the relationships of the evoked hemodynamic and electrophysiological responses particularly with regards to nonlinarites in the neural response for interstimulus intervals between 50 to 500 ms. Specifically, the magnitude of the evoked neural and hemodynamic signals in a region-of-interest analysis for the parametric pulsed-pair response was compared.

Subject Population
In this study, a total of 22 subjects (15 males, 7 females; age range 21 to 45 years; all right handed) were recruited for functional imaging studies. All subjects provided written informed consent prior to study procedures. This study was approved by the University of Pittsburgh Institutional Review Board.
The study consisted of two separate scanning sessions for the MEG-fNIRS and MRI-fNIRS. The MEG-fNIRS scan was done first in all cases due to the potential loss of subjects whose head size or body dimensions precluded them from MEG scanning. A total of four subjects was excluded prior to scanning because their head would not fit inside the MEG scanner while wearing the fNIRS sensors. A further two subjects failed to yield adequate fNIRS signals after initial setup and were excluded. The remaining 16 subjects were scanned with concurrent fNIRS-MEG. Of the 16 remaining subjects, 11 were additionally scanned inside the MRI with concurrent fNIRS. Five subjects were lost to follow-up or withdrew from the study after the first visit. Of the 11 subjects that completed both the MEG and fMRI sessions, a further 4 were rejected due to excessive eye-blink or head motion artifacts in the MEG analysis. Thus, a final total of 7 subjects (4 males, 3 females) were used in the final analysis of the study having completed all aspects of the protocol.

Experimental Protocol
Unilateral median nerve stimulation was conducted using a computer controlled GRASS S88X stimulator (Natus Neurology Inc., Warwick, Rhode Island). A stimulating bar electrode was positioned on the right median nerve and the voltages were adjusted to match the motor threshold as judged by twitching of the thumb on the stimulated hand. The stimulating voltage was then set at the determined motor-threshold for the study. The median nerve stimulation had duration of 20 μs with a monophasic pulse. Two experimental tasks were performed. First, a functional localizer task consisting of repeated 10-s blocks of stimulation at 4 Hz with a 10-s interblock rest period (5 min total time; 15 blocks) was administered. The purpose of this task was to provide spatial localization information to define regions-of-interest to the pulsed-pair stimulation task. Second, several 5-min scans of event-related pulsed-pair stimulation were then given. The pulsed-pair condition consisted of a pair of two median nerve pulses separated by an interstimulus interval of 50, 100, 150, 200, 300, 400, or 500 ms. The stimulus conditions were presented in randomized order with a minimum spacing of 8-s between events to avoid potential nonlinearity in the hemodynamic response based on the work by Cannestra et al., 16 which suggested a minimum of 4-s interval was required. A total of 32 pairs of stimulus events was delivered in each 5-min scan. This scan was repeated 4 to 6 times per subject depending on subject tolerance and time constraints. The experimental procedure for both the MEG-NIRS and MRI-NIRS portions of the study were identical. An MEG-compatible bar electrode obtained from Elekta Neuromag (Helsinki, Finland) was used in this study. An identical duplicate electrode was used for the MRI portion of the study, which was verified to have negligible RF-heating and no ferromagnetic characteristics.
In order to control the stimulus timing, a custom dynamically linked library (DLL) was written to digitally interface to the GRASS stimulator via a universal serial bus cable. This allowed digital control of the frequency and duration of the stimulation with microsecond precision. An experimental paradigm was written in Eprime-2 (Psychology Software Tools; Pittsburgh, Pennsylvania) to control the GRASS device via the custom DLL program to prime the stimulator with the correct pulse interval several seconds prior to stimulation. The start of each stimulus event was triggered by Eprime and the GRASS device then controlled the timing of the second pulse. In this way, the pulsed-pair stimulus was delivered with better than millisecond precision as measured by the MEG auxiliary channel. The stimulus output from the GRASS device was sent to both the fNIRS and MEG systems for recording. A current isolation device was used to limit stimulation of the subject.

Analysis Methods
A schematic overview of our complete multimodal analysis pipeline is shown in Fig. 1. A key feature of this pipeline is the use of the subject anatomy and cortical surface model generated by the FreeSurfer tools available from the Massachusetts General Hospital. 17 The extracted cortical surfaces, which are registered among subjects to allow group analysis, are used in all three functional imaging modalities, which allows direct comparisons of each modality. The individual components of this pipeline are described as follows.

MEG Acquisition and Analysis
MEG data were collected using a whole-head 306-channel Elekta Neuromag ® Vectorview system (Elekta Neuromag Oy, Helsinki, Finland). Concurrent electrooculogram (EOG) was also acquired. This system is composed of 102 triple-sensors consisting of a magnetometer, longitudinal gradiometer, and a latitudinal gradiometer and is contained inside a magnetically shielded room (MSR). During the acquisition, subjects were inside the MSR (∼4.6 m × 3.6 m × 3 m in dimension) sitting comfortably with their head positioned in a helmet containing the sensors (Fig. 2). MEG scans were sampled at 1000 Hz and subsequently down-sampled offline to 250 Hz. The MEG channels were manually inspected and obviously bad channels were removed from further analysis. The MEG signals were then bandpass filtered using cutoffs of 1 and 50 Hz to remove the 60 Hz line noise and slow drifts. Temporal spatial signal separation (tSSS) was employed to remove external magnetic field contributions in the signal recordings, which can lead to spurious signals that appear to have originated from physiological changes but are in fact artifacts external to the MEG sensors. 18 Data were then averaged into trials and trials with excess eye blinks or eye movements, as assessed by the maximal difference in the EOG over the trial, were rejected. Data from each condition were subsequently averaged offline.
For source localization, the cortically constrained minimum norm estimate (MNE) inverse solution was employed using the MNE ™ software suite v2.7. 19 Briefly, a linear inverse operator W is applied to the measured signal to calculate the MNE E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 4 ; 3 2 6 ; 4 3 7 yðtÞ ¼ WxðtÞ; where xðtÞ represents the MEG channel data at time t and yðtÞ is the corresponding current projected onto the cortical surface. The expression of W is calculated using the L 2 norm, which yields where A is the free source orientation solution of the forward problem calculated using the boundary element method. C and R are the noise and source covariance matrices, respectively. The noise covariance matrix was calculated using intertrial time. A defines the transformation from orthogonal unit current dipoles to measured magnetic fields and uses the boundary element method. 20,21 λ is a weighting factor that is used to avoid the magnification of errors in the data and λ 2 ≈ ð1∕SNRÞ. We used a value of 3 for this, as is typically done in MEG analysis (MNE software user's guide version 2.7, 2009; 22 ). Furthermore, because cortical neurons are known to be preferentially oriented perpendicular to the cortical surface, we used a loose orientation constraint. Specifically, the components transverse to the surface were multiplied by 0.4. To compensate for the bias toward superficial currents of the MNE, a scaling factor (i.e., depth weighting) of 0.8 is applied to R. 23 For group-level maps, the data from each participant was morphed onto a common brain. Statistical values were assessed using a paired t-statistic across the seven participants for the MEG response to the median nerve stimulation versus prestimulation baseline. The peak time was chosen individually for each subject based on the time that showed the greatest change relative to baseline averaged across all MEG sensors (mean ¼ 74.4 ms; std ¼ 15.5 ms). For the ROI analysis, ROIs were drawn at the individual subject level. The ROI corresponding to primary somatosensory cortex was determined from the response in the functional localizer scan relative to prestimulus time and applied to the pulsed-pair condition.

MRI Acquisition and Analysis
Structural and functional MRI was performed using a 3T Siemens TIM TRIO scanner using a 12-channel head coil. A structural T1-weighted image was acquired (Siemens' MPRAGE; 176 × 256 × 192 1 mm isotropic; TE ¼ 3.52 ms; TR ¼ 2300 ms). A total of 4 to 6 5 min BOLD EPI scans were performed (64 × 64 × 31 3.5 mm isotropic; TR ¼ 2000 ms; TE ¼ 30 ms; FA ¼ 90 deg). Functional MRI analysis was performed using the AFNI and SUMA toolset from the NIH. 24,25 In brief, functional EPI data were coregistered to the structural volume and motion corrected to the third TR. The slices were temporally aligned and volumes were blurred using the default 4-mm FWHM kernel. General linear model analysis based on the timing of the seven pulsed-pair and localizerconditions was performed using the default gamma-variant canonical basis function (GAM) 1 with additional motion regressor terms. Group level ANOVA was performed within AFNI.
The structural T1-weighted anatomical MRI image for each subject was processed using the FreeSurfer pipeline in order to generate the segmented anatomical models and registered cortical surface models used for fNIRS, MEG, and fMRI analysis procedures. Following registration and brain segmentation using Freesurfer, the MNE toolkit 22 was used to additionally create boundary element models for the skin and outer/inner skull, which allowed construction of a 5-layered head model (skin, skull, CSF, gray matter, and white matter) for each subject which was used to generate the fNIRS forward model using the optical properties described in Fang and Boas 26 for the 690-and 830-nm wavelengths. The Monte Carlo program tMCimg 27 was used to simulate the optical forward model. The same boundary element and brain surface models were used to generate the MEG forward model as part of the MNE MEG analysis pipeline developed by Hamalainen et al. at the Massachusetts General Hospital. 22 The extracted white matter and pial brain surfaces from FreeSurfer were used in the inverse models for both fNIRS and MEG and for the SUMA-AFNI 24,25 surface analysis of the fMRI data. In this way, the estimates of brain activation from all three modalities were derived in the same cortical surface space for direct comparisons of the results.

fNIRS Acquisition
For both MRI and MEG portions of the study, fNIRS signals were recorded using a commercial TechEn CW6 (Milford, Massachusetts) continuous-wave fNIRS system. This system contains up to 16 laser diodes each at 690 and 830 nm (12 and 8 mW, respectively), which are coupled through fiber optic cables mounted onto a head cap. The head cap was constructed of Velcro and plastic materials and was custom built for this study. Due to the limitations of head-clearance inside the MEG scanner, the fNIRS cap had to be restricted to eight detectors and four source positions containing both wavelengths on only one hemisphere of the head. The same design was used for both parts of the study although dedicated identical sets of fiber optics and head caps were used to avoid magnetic contamination of the MEG. The probe extended from left frontal cortex (near 10 to 20 coordinate "F3") to posterior partial (near 10 to 20 coordinate "P3") with a source-detector spacing of 3.2 cm as shown in Fig. 2.

MEG-fNIRS Setup
The fNIRS cap used in this study had a profile height of 11 mm, which reduced the spacing between the MEG SQUID detectors and the subject's head. Four head position indicator (HPI) coils were glued to the subject's scalp using collodion gel. These were used to coregister the fNIRS and MEG sensor positions after positioning the fNIRS cap on the subject's head. Two electrodes were positioned above and below the left eye and above the corner of the left eye and below the corner of the right eye to monitor eye blinks and saccades. Cardiac activity was recorded on an ECG channel using two electrodes that were attached just below the collarbone of the subjects. Prior to MEG scanning, a Polhemus Isotrack (Polhemus, Inc., Colchester, Vermont) positioning system (integrated with the MEG acquisition software) was used to register the HPI coils and to landmark fiducial points (nasion and preauricular points) on the subject's head. These were then used to register both the fNIRS and MEG sensor locations to the subject's anatomical MRI data. Based on previous work, we expect the localization error of this method to be around 3 to 5 mm. 28 Additional points (up to 100) were digitized on the subject's head and face to help facilitate coregistration of MEG data with the structural MRI's. Immediately prior to MEG acquisitions, the HPI coils were energized and their positions localized through the MEG acquisition interface. This registration information was then extracted from the MEG data files and combined with the location of the fNIRS sensors from the isotrack data using a rigid-body transform to localize the fNIRS sensors on the registered MRI volumes from the second session.
For fNIRS recordings, an MEG-dedicated set of fiber optics (10-m long) was passed through the access port (wave-guide) in the side of the MSR. The fNIRS instrument was located next to the MEG acquisition computers. The stimulus timing from the median nerve device was sent to both the MEG and fNIRS instruments through auxiliary recording channels built into both systems. Because the fNIRS instrument was sampled at a lower rate, a 555-timer pulse-elongating circuit had to be constructed to elongate the 20-μs stimulating pulse to 500 ms for recording by the fNIRS system.

MRI-fNIRS Setup
Functional NIRS collection during the MRI portion of the study was similar to the MEG collection with the exception that a MRI-dedicated set of 10-m fiber optics and head cap was used to avoid magnetic contamination to the MEG system. The CW6 instrument was positioned in the control room for the MRI scanner and the fiber optics were run through the MR waveguide into the scanner room. The fNIRS head cap used in the MRI portion of study was identical to the one used in the MEG portion of the study. The head cap contained vitamin E capsules in place of the four MEG HPI coil positions, which were used for registration of the cap from the T1weighted MRI image for each subject. In previous work, we have estimated the registration error of this approach to be about 5 mm. 14 The timing of the stimulus paradigm was triggered from the MRI scanner.

fNIRS Analysis
Analysis of fNIRS results was performed using custom MATLAB scripts. Raw 20-Hz fNIRS data were converted to optical density (absorption) at both 690 and 830 nm for all measurement pairs. The unfiltered optical density data for all measurement pairs were then analyzed using a statistical general linear model as described in Barker et al. 11 This model uses an iterative estimation of an autoregressive (AR) whitening filter and weighted-least-squares regression to reduce the effects of serially correlated noise and motion artifacts in the fNIRS data. In brief, a regression model (design matrix; X) was constructed from convolution of a canonical hemodynamic response function (default AFNI GAM 1 function; 24 ) with the stimulus onsets and durations of the median nerve stimulation events. Thus, a linear regression model is constructed of the form Y ¼ X · β, where Y is the vector of time-points for each optical measurement pair (optical density) and β is the estimated weights of the regressors. The values of β are first estimated from a weighted least-squares procedure, which reduces the influence of noise outliers (robustfit function in MATLAB 29 ).
A pth-order AR model is then fit to the residual error of the regression fit where p is estimated by an AIC criteria search. The AR coefficients are used to construct a whitening filter (S −1 ), which is applied to both the left and right hand sides of the regression model (e.g., S −1 · Y ¼ S −1 · X · β). The regression model is again solved using robust regression methods and the procedure is iterated until convergence (usually within 2 to 3 iterations). This is repeated for each optical measurement channel and wavelength. Our previous work in developing this model, described in Barker et al., 11 had shown that this model was robust to both spike and shift types of motion artifacts in the fNIRS data. In particular, the presence of infrequent, large amplitude motion artifacts in the fNIRS data results in a heavy-tailed statistical distribution of noise, which can be addressed with robust regression methods to treat these outliers. In addition, we had shown that correction for serially correlated errors due to systemic physiology and drifts produced much more accurate estimates of model errors and type-I error rates compared to least squares estimation. Thus, this model is fairly robust to false-discoveries due to superficial physiological contamination. 11 This model performs best when raw (unfiltered) data are used since removing portions of the frequency spectrum can make the prewhitening procedure ill-posed and the algorithm numerically unstable. The outputs of the regression models (e.g., estimates of β and the estimate of the covariance errors in β associated with the optical density changes for each fNIRS measurement channel) were then reconstructed into estimates of brain activity. For image reconstruction, we used a random-effects cortical surface topography model, which had been previously described in Abdelnour et al. 13 In this model, this fNIRS forward model (sensitivity model) is computed on the cortical surface of the brain using the extracted surfaces from the Freesurfer model described in the previous section. Of note, this same cortical surface is also used in the MEG and fMRI models via AFNI-SUMA. The registered "spherical" inflated surface for each hemisphere from the Freesurfer model is reparameterized using spherical wavelet basis set described in Abdelnour et al. 14 Spherical wavelets are used to describe the image (e.g., brain activity map) on the surface of the cortex at different spatial resolutions similar to more conventional wavelet representations with the distinction that these spherical wavelets are specifically designed to model a two-dimensional topography on a closed three-dimensional surface. The spherical wavelets are derived from a recursive subdivision of an icosahedral mesh, which can easily be extracted from the existing Freesurfer cortical surfaces from the icosahedral levels 1 to 4 (2562 nodes per hemisphere with about 6.2-mm vertex spacing and 39-mm 2 face areas). Since the registered "spherical" surfaces are used in the FreeSurfer model for each subject, the spherical wavelets map on to approximately the same anatomical regions across subjects. For example, the nth parameter in the model represents brain activation at the same cortical gyrus for all subjects even though that gyrus may be in a slightly different Cartesian location in the space of the pial brain surface. Of further note, this surface mapping feature of the Freesurfer models allows activations to be transferred and averaged among subjects and is used in both the AFNI-SUMA and MEG-MNE pipelines for group level models. We note that this is the same convention and mesh used in the MEG MNE reconstruction through the Freesurfer pipeline.
Because the spherical wavelet model parameterizes each individual subject's fNIRS forward model into the same registered brain space, the models for all subjects can be concatenated into a single large random-effects inverse model and solved simultaneously for all subjects. Our random-effects group level fNIRS model was previously described in Abdelnour et al. 13 In brief, this group level analysis model is similar to the group-level image reconstruction of MEG data described in Mattout el al. 30 The fNIRS forward model for all individual subjects (L i ) can be concatenated into a singleinverse operator such that    In this random-effects model, each individual subject's image is modeled as a perturbation on the group average (e.g., β 1 ¼ β Group þ Δβ 1 ). Because the fNIRS inverse model is ill-posed (meaning that there are a large number of equally plausible solutions that model the data), this random-effects approach attempts to find an estimate of the group average that is self-consistent with all the individual subjects. Our previous work 13 had shown that this approach resulted in a much smaller point-spread function and higher effect sizes in group-level analysis compared to a more conventional twostep approach of reconstructing each individual subject and then averaging the results. While this approach has also been shown to be effective in MEG analysis, 30 this is not yet implemented in the MEG pipeline used by MNE and thus, it is important to point out that our MEG processing still used this two-step approach.
The fNIRS random-effects inverse model was solved using a hierarchical Bayesian model as described by Abdelnour et al. 13 Restricted maximum likelihood and expectation-maximization recursion was used to estimate the noise hyper-parameters of the model. The covariance of the regression coefficients from the general linear model (β) was used to construct the covariance prior on the measurement noise, whereas a scalar identity operator was used to model oxy-and deoxy-hemoglobin priors. Thus, a total of three hyper-parameters (λ) are estimated in this model. A description of the algorithm is found in Abdelnour et al. 13 and follows the description of the hierarchal inverse model used in the MEG and EEG image analysis in SPM. 30,31

Results
The fMRI, fNIRS, and MEG modalities all showed significant changes in brain activity during the median nerve localizer task. Figure 3 shows the group level results for the three modalities. In all modalities, the activation was located in the left somatosensory cortex as shown in Fig. 3(a). The fNIRS results show the group-level results from across both imaging sessions. A closer view of the somatosensory cortex is shown in Fig. 3(b), which shows the location of all three modalities along the primary somatosensory gyrus. Figure 4 shows the results of the fNIRS image reconstructions for the fMRI-fNIRS and MEG-fNIRS sessions independently to examine the intersession differences in the fNIRS signals. The fMRI-fNIRS and MEG-fNIRS sessions show similar regions of the left somatosensory region particularly for the oxy-hemoglobin images. The MEG-fNIRS session data showed a bit more noise for the deoxy-hemoglobin image compared to the fMRI-fNIRS session or combined data likely related to lower data quality due to compression of the fNIRS cap within the MEG sensor helmet.
Based on the image reconstructions for the three modalities, the center-of-mass of the activation patterns for all images are given in Table 1. The location of the center-of-mass of each modality is shown in Fig. 5(a). The displacement between the centers ( Table 2) of the fMRI-BOLD signal was 26, 20, and 19 mm from the MEG, fNIRS oxy-hemoglobin, and fNIRS deoxy-hemoglobin cluster centers, respectively, when examined over the entire brain. This result was consistent with the recent fMRI-MEG study by Stevenson et al. using similar median nerve tasks that showed a displacement between the dipole MEG model and the BOLD signal of 19 AE 6 mm. 32 The N20 component was used to localize the MEG response from the localizer median nerve task. The center of the MEG was 27 and 29 mm from the oxy-and deoxy-hemoglobin centers. Finally, the center of the oxy-hemoglobin cluster was 3 mm from the center of the deoxy-hemoglobin cluster. Due to both MEG and fNIRS having greater sensitivity to shallower sources, we hypothesized that much of the relative displacement across modalities was due to discrepancies for deeper sources. To test this, we examined the center of mass displacement as a factor of whether or not deeper sources were included in the map [ Fig. 5(b)]. This analysis demonstrated that there is strong correspondence across the modalities at shallower depths, peaking at 9 mm (HbO 2 ), 16 mm (Hb), 10 mm (MEG) from the fMRI-BOLD center of mass when including the image up to 20 mm from the scalp. As deeper sources were included in Fig. 3 Reconstructed images of brain activation to the median nerve localizer task for the fMRI, MEG, and fNIRS modalities. The T -score for the localizer median nerve task is shown for each modality and thresh-held at p < 0.05. (a) The reconstructed activation maps for the fMRI-BOLD, MEG, and fNIRS signals (average of both the fMRI and MEG sessions) are shown. Activation was observed on the left somatosensory (S1) region of the brain. the analysis, the displacement of the centers of mass converges to what is shown in Table 2.
In Fig. 5(c), we show the depth of the activation from the surface of the scalp for each imaging modality. The massweighted average depth of the activation from the nearest point on the surface of the scalp as detected by fMRI-BOLD was 35 mm. For MEG, this depth was about 7 mm more shallow at 28 mm and for fNIRS this was almost 15 mm shallower at 19 mm from the scalp. Fig. 5(c) shows the histogram plots of the depth location of the detected activation areas (p < 0.05) for each modality. For fNIRS, ∼30% to 35% of the reconstructed signal was between 18 to 22 mm from the surface of the scalp. There were no signals reconstructed beyond about 28 mm.
In contrast, for the fMRI-BOLD, roughly 50% of the signal was located at a depth of 28 mm or greater, which means that the fNIRS method missed about half of the signal detected by BOLD because of the limited fNIRS depth penetration. The depth of the MEG reconstructed signal was deeper than the fNIRS, but still had a fall off around 30 mm compared to the BOLD signal.
To further quantify the relationship among the modalities, we examined the spatial correlation between each modality pair as a factor of whether or not deeper sources were included in the correlation. In Fig. 5(d), we show the correlation of the fNIRS and MEG signals as a function of the inclusion distance from the surface of the scalp. The correlation for both MEG and The T -score for the localizer median nerve task is shown. fNIRS with the BOLD signal was strongest only if the superficial cortical points were included. As the inclusion depth increases, the correlations decrease to approach the wholebrain asymptotic limit, which includes comparisons in deep subcortical areas where neither MEG nor fNIRS have strong sensitivity. Both the whole-brain and peak correlation at depth are provided in Table 3. The whole-brain limit of correlation between fMRI-BOLD and the MEG was R ¼ 0.17. Likewise, Fig. 5 The center-of-mass of the localizer activation patterns was calculated for the four reconstructed image types and two NIRS sessions. There was a slight displacement of the fNIRS and MEG signals toward the cortical surface compared to the center of BOLD activation. In addition, the fNIRS results from the fNIRS-MEG session were slightly more frontal compared to the other data types. The data are additionally summarized in Tables 1 and 2. Table 2 Center-of-mass activation displacements. The displacement in millimeters between the center-of-mass of the estimated localizer task responses for each modality is shown. The upper-right half of the table (italics) shows the minimum displacement between each modality along the depth from the scalp surface [e.g., also see Fig. 5(b)]. The bottom-left half (not italics) shows the center-of-mass displacement over the entire brain, including the subcortical regions. Distances cited in millimeters. the correlation between BOLD and oxy-and deoxy-hemoglobin was R ¼ 0.28 and R ¼ −0.25, respectively. The correlation between the MEG and oxy-and deoxy-hemoglobin was R ¼ 0.16 and R ¼ −0.14, respectively. However, by restricting the correlation analysis to the superficial cortex only, the maximum correlation between the MEG and fMRI-BOLD was found to be R ¼ 0.54 and was R ¼ 0.57 and R ¼ −0.48 between the BOLD and fNIRS oxy-and deoxy-hemoglobin signals, respectively (Fig. 5). Likewise, the maximum correlation between the MEG and fNIRS reconstructions was 0.46 and 0.37 for oxy-and deoxy-hemoglobin. The maximum correlation between fNIRS oxy-and deoxy-hemoglobin was R ¼ −0.94. The spatial correlation results between all other combinations are provided in Table 3. All correlations (including when deep sources are included in the analysis) were statistically significant (p < 0.001).

Pulsed Pair Results
In the second half of the study, we performed an event-related pulsed-pair median nerve study. The pulsed-pair task uses two stimuli placed 50 to 500 ms apart. A slow event-related design with a minimum 8-s ISI was used to avoid additional nonlinearity in the hemodynamic response unrelated to the neuralrefractory effect. 16 As shown in Fig. 6, when the second pulse was <200 ms after the first, the MEG ERP response was significantly diminished by the neural refractory period. At 300 ms and later, the second ERP response was again visible. A similar task had previously been used by Ogawa et al. 33 showing a similar effect in rodent forepaw stimulation and was more recently examined for multimodal fMRI-MEG by Stevenson et al. 34 Using the region-of-interest defined as the union of activations from all of the modalities observed in the localizer task the magnitude of brain activity for all seven interpair intervals was calculated from each modality. In Fig. 7, grand average of the normalized magnitude of the MEG [area-under-curve ðAUCÞ ¼ P 600 ms t¼0 jERPðtÞj], fMRI, and fNIRS signals for the pulsed-pair conditions is presented. The signals are shown normalized to the mean of the seven conditions. We found significant correlations of R ¼ 0.80, 0.52, and −0.70 between the magnitude of the fMRI and fNIRS [oxy-/deoxy-hemoglobin] responses compared with the AUC of the MEG response using colocalized regions-of-interest. The BOLD and fNIRS magnitudes were also significantly correlated for oxy-and deoxy-hemoglobin responses (R ¼ 0.27 and −0.89, respectively).
In addition, the amplitude of the evoked pulsed-pair response was also estimated using modality-specific regions-of-interest. For the fNIRS data, the region-of-interest was defined in fNIRS channel (source-detector) space as the average value over all statistically significant channels (p < 0.05) on the probe. Using the modality-specific regions-of-interest, we still found statistically significant (p < 0.05) correlations between the MEG signal (AUC) and the fNIRS channel oxy-and deoxyhemoglobin data from the MEG-fNIRS session (R ¼ 0.75 and −0.69, respectively). The fNIRS channel-based correlations to the MRI-BOLD signal went for down to R ¼ 0.26 and −0.46 for oxy-and deoxy-hemoglobin. These lower correlations observed in fNIRS channel space compared to the reconstructed images may be due to errors introduced from intersubject registration and anatomy, which are accounted for in the image reconstruction process but not in the raw channel space data.

Discussion
In this work, we examined the spatial localization of reconstructed fNIRS signals and compared this to both fMRI and MEG measurements. Although previous work has established the general correspondence of these measurements from a physiological and biophysical basis (reviewed in Ref. 10), it is unclear how well spatial agreements can be achieved in practice given the inherent difficulties and ill-posed nature of the reconstruction of low-density fNIRS measurements into spatial images. Previously, we have described and numerically validated several advanced methods for fNIRS image reconstruction that incorporate anatomical information from structural MRI, 14 group-level random effects models, 12 and hierarchical regularization models. 13 In addition, we have recently detailed several methods for improved time-series analysis and generalizations of the linear model to deal with fNIRS specific noise structures. 11,35,36 In this work, we have presented an analysis pipeline combining these recent advancements that allows all Table 3 Spatial correlation. The correlation of the activation patterns between the estimated localizer responses. The upper-right half of the table (italics) shows the maximum correlation values between each modality along the depth from the scalp surface (e.g., also see Fig. 5). The bottom-left half (italics) shows the correlation over the entire brain, including the subcortical regions that showed activity in fMRI but not fNIRS nor MEG. three modalities (MEG, fMRI, and fNIRS) to be presented and compared in the same cortical brain space.

Proposed Analysis Pipeline
As detailed in Fig. 1, our analysis of all three modalities was centered upon the use of the inflated cortical surface from FreeSurfer 17 as the registration and reconstruction basis for each modality. In the case of both fNIRS and MEG, this cortical surface was used to restrict the image reconstruction problem to a two-dimensional (manifold) estimation problem directly on the folded three-dimensional cortical surface. Similarly, this FreeSurfer cortical representation is used in the MNE 22 package for MEG analysis and as part of the AFNI/SUMA analysis package. 24,25 The use a common basis upon the cortical surface for analysis in all three modalities allowed direct comparisons of the estimated brain activation patterns across modalities. However, we also propose that this approach can be useful in future work to better utilize the mutual information among multiple modalities. For example, in our previous work in Adbelnour et al., 13 we described how maximum likelihood methods could be used to introduce multiple regularization priors for fNIRS image reconstruction based on prior fMRI (or MEG) data. In this approach, multimodal data are used to constrain the reconstruction of the fNIRS data as a statistical  prior. Since all three modalities are defined in the same spatial basis support, this pipeline will allow future use of such methods.

Comparision of fNIRS, fMRI, and MEG Signals
The spatial reconstruction of fNIRS signals requires a solution to the optical inverse problem (e.g., estimation of brain signals from a limited number of surface measurements). Therefore, spatial comparisons of fNIRS and other modalities are sensitive to the means of solving this inverse problem. To get around this issue, previous work by Sassaroli et al. 37 and Huppert et al. 38 both looked at the consistency of the measured sensor-space fNIRS data with the fMRI image by using the fNIRS sensitivity model to forward-project the fMRI into the measurement space of the fNIRS signals for comparison and thereby avoiding the mathematical image reconstruction process. These works showed that the fNIRS channel-data were spatially consistent with the fMRI data, but did not address the question of how well fNIRS on its own can localize signals within the brain given the current state-of-the-art in data collection and image reconstruction methods for fNIRS imaging. Using high-density optical tomography, Eggebrecht et al. 39 showed a spatial correlation between NIRS and fMRI of between R ¼ 0.26 (verb generation task) and R ¼ 0.83 (reading task) at similar 6-mm fMRI smoothing as used in our current study. However, in contrast to our study, Eggebrecht et al. used a high-density NIRS system with a total of up to 644 NIRS source-detector pairs in comparison to the standard-density probe with 28 pairs used in this current study. Given the lower cost, ubiquity, and increased flexibility of standard density acquisition as was used in this study, it is notable that the correlations reported here compare favorably with those reported in the high-density fNIRS study by Eggebrecht et al. mentioned above. This is particularly relevant to multimodal imaging, as high-density NIRS systems are not compatible with concurrent data collection. More recently, the work by Cui et al. 40 looked at both sensor-space and imagebased comparisons of fNIRS and fMRI in a variety of cognitive and sensory-motor tasks and found spatial correlations of R ¼ 0.23 −0.26 based on a simple spline interpretation to display the fNIRS data into a standardized brain template. The relationship between fNIRS and electrophysiology recordings has also been examined. The review by Shibasaki 41 discusses the topic of correlation between the hemodynamic and electrophysiological signals, including a discussion of fNIRS. The localization of functional tasks in concurrent MEG-fNIRS and EEG-fNIRS has also been examined by a few groups. Sander et al. 42 looked at direct current (DC) MEG and fNIRS signals during a motor task study. A similar study was also done by Mackert et al. 43 that also showed qualitative agreement between DC-MEG and fNIRS. Similarly, Seki et al. 44 looked at auditory function using a custom built unshielded MEG and fNIRS system. A more comprehensive study using conventional MEG was done by Ou et al. 45 Ou et al. used a median nerve electrical stimulation task to look at the relationships of the fNIRS and MEG signals under a parametric stimulation task. They examined the relationship of the amplitude of the optically measured oxy-and deoxy-hemoglobin signals with respect to various physiological peaks in the MEG signal. Specifically, they found that the amplitude of the N20 and P60 peaks best predicted the change in the hemodynamic signal for changes in the duration of the stimulation task. This study also looked at the spatial localization of the fNIRS and MEG signals and found spatial agreement with MEG source localization using a dipole model and the nearest optical channels (e.g., source localization was not used in fNIRS).
Using our proposed pipeline, we found close agreement with all three modalities as presented in Fig. 3. The activations for all three modalities were localized to the primary somatosensory (postcentral) gyrus and primary motor (precentral) gyrus. Although the fNIRS and MEG spatial resolution is generally coarser than this, the additional information added by the cortical surface used in the image reconstruction process adds spatial information to constrain the images. Thus, with the addition of this structural information, we found that the overall spatial resolution of the fNIRS images was comparable to the BOLD fMRI. We found that for shallower regions of the brain, there was only a 9-and 16-mm displacement between the BOLD and fNIRS activations [center-of-mass; HbO 2 and Hb, respectively)] and R ¼ 0.47 and R ¼ −0.42 correlation between the two maps at a depth of 20 mm. This error was similar to that observed between the MEG and BOLD signal (10.4 mm) with a similar correlation (R ¼ 0.43). To our knowledge, this is the first study to use the pulsed-pair task with fNIRS. However, our results are similar to the work by Stevenson et al. 34 that used a similar pulsed-pair median nerve task during independent sessions of fMRI and MEG in seven subjects. They reported good spatial agreement between the two modalities. In their study, the displacement of the average BOLD and MEG (dipole) locations was found to be 19 AE 6 mm, which is similar to the 25-mm displacement found in this study for the whole head analysis, but greater than the 7.1-mm displacement found when only the superficial layers of the cortex were considered in the analysis. The relatively large inaccuracy of the signals in MEG and fNIRS for depths beyond 30 mm suggests that it is prohibitively difficult to image deep brain sources, such as subcortical regions, with these methods. Previously, this falloff of signal with depth had been modeled using simulated data, [46][47][48] but the current work provides an empirical examination of this issue. One caveat of the current study is that it is done in the context of median nerve stimulation, which provides a relative context for comparing localization.

Test-Retest Reliability
Since the fNIRS data were recorded in two separate sessions (once with fMRI and once with MEG), we also examined the test-retest reliability of the fNIRS signals and image reconstruction. As was presented in Fig. 4, we found close agreement across the sessions for the oxy-hemoglobin estimates of the response to the localizer task. The deoxy-hemoglobin images were slightly noisier, particularly in the case of the MEG-NIRS session data, but were still strongly correlated with the fMRI-NIRS session data (Table 3). This indicates that the fNIRS results were spatially consistent between the two sessions. There was a slight frontal bias of the estimated activation in the MEG-NIRS session in which more activation was observed in the frontal lateral region of the probe (see Figs. 4 and 5).

Conclusions
The results of this study show close agreement in spatial localization of fNIRS with MEG and fMRI results. In addition, the magnitude of estimate brain responses for these three modalities was linearly correlated during a parametric task. The main source of difference across the modalities was the decreased sensitivity with respect to depth seen in fNIRS and MEG. These results validate the concurrent recording approach for future studies of neurovascular coupling and multimodal fusion with fNIRS building a bridge between fMRI and MEG.

Disclosures
None of the authors has any conflicts of interest to disclose.