Impact of Anatomical Variability on Sensitivity Profile in fNIRS–MRI Integration

Functional near-infrared spectroscopy (fNIRS) is an important non-invasive technique used to monitor cortical activity. However, a varying sensitivity of surface channels vs. cortical structures may suggest integrating the fNIRS with the subject-specific anatomy (SSA) obtained from routine MRI. Actual processing tools permit the computation of the SSA forward problem (i.e., cortex to channel sensitivity) and next, a regularized solution of the inverse problem to map the fNIRS signals onto the cortex. The focus of this study is on the analysis of the forward problem to quantify the effect of inter-subject variability. Thirteen young adults (six males, seven females, age 29.3 ± 4.3) underwent both an MRI scan and a motor grasping task with a continuous wave fNIRS system of 102 measurement channels with optodes placed according to a 10/5 system. The fNIRS sensitivity profile was estimated using Monte Carlo simulations on each SSA and on three major atlases (i.e., Colin27, ICBM152 and FSAverage) for comparison. In each SSA, the average sensitivity curves were obtained by aligning the 102 channels and segmenting them by depth quartiles. The first quartile (depth < 11.8 (0.7) mm, median (IQR)) covered 0.391 (0.087)% of the total sensitivity profile, while the second one (depth < 13.6 (0.7) mm) covered 0.292 (0.009)%, hence indicating that about 70% of the signal was from the gyri. The sensitivity bell-shape was broad in the source–detector direction (20.953 (5.379) mm FWHM, first depth quartile) and steeper in the transversal one (6.082 (2.086) mm). The sensitivity of channels vs. different cortical areas based on SSA were analyzed finding high dispersions among subjects and large differences with atlas-based evaluations. Moreover, the inverse cortical mapping for the grasping task showed differences between SSA and atlas based solutions. In conclusion, integration with MRI SSA can significantly improve fNIRS interpretation.


Introduction
Functional near-infrared spectroscopy (fNIRS) is gaining relevance in the imaging of cortical activity in a variety of tasks [1][2][3][4]. Compared to functional MRI (fMRI), it allows for the performance of low-cost and ecological measurements of brain activity in an open environment, also minimizing discomfort to participants. Furthermore, fNIRS provides a direct estimation of both the concentration changes in oxygenated HbO 2 and deoxygenated hemoglobin (HbR), hence providing deeper insight to the hemodynamic response to neural activation compared to the blood oxygen level-dependent signal of fMRI. Moreover, fNIRS is well-suited to multimodal integration with other neuroimaging techniques such as EEG [5] and MRI [6,7].
Acquisitions are based on a grid of optodes alternating near-infrared light sources and detectors, whose neighboring pairs form a measurement channel. However, fNIRS' practical advantages and noninvasiveness requirements for clinical applications must face two major technological issues: (i) the need of accurate signal processing to reduce

Dataset and Experimental Set-Up
At IRCCS Fondazione Don Carlo Gnocchi, 13 healthy young adults (6 males, 7 females, age 29.3 ± 4.3) underwent both fNIRS acquisitions and anatomical MRI scans. The experiment was approved by the IRCCS Fondazione Don Carlo Gnocchi ethical committee, and all volunteers provided their informed consent. An expert neuroradiologist examined each MRI to exclude the presence of any pathology.
The focus of this study is on the analysis of the forward problem and the impact of SSA variability, ahead of the regularization implied by the inverse solution. Nonethe-Sensors 2023, 23, 2089 3 of 19 less, application examples of fNIRS cortical mapping are shown on representative subjects performing a motor task. Data were acquired using a continuous-wave system at λ = 760 nm and λ = 850 nm wavelengths (NIRScoutX 32 × 32, NIRx Medizintechnik, Berlin, Germany) that employs 32 LED sources and 32 avalanche photodiode detectors placed according to the international 10/5 EEG system (Figure 1). Source and detector pairs were combined into a total of 102 measurement channels.
The focus of this study is on the analysis of the forward problem and the impact of SSA variability, ahead of the regularization implied by the inverse solution. Nonetheless, application examples of fNIRS cortical mapping are shown on representative subjects performing a motor task. Data were acquired using a continuous-wave system at λ 760 nm and λ 850 nm wavelengths (NIRScoutX 32 × 32, NIRx Medizintechnik, Berlin, Germany) that employs 32 LED sources and 32 avalanche photodiode detectors placed according to the international 10/5 EEG system (Figure 1). Source and detector pairs were combined into a total of 102 measurement channels. The functional task consisted of a block-designed motor-grasping task, where participants had to alternatively move their left or right hand [8]. The fNIRS signals were preprocessed using the Nirstorm package (https://github.com/Nirstorm/nirstorm (accessed on 5 February 2023)) in Brainstorm [34] to extract the relative changes in oxygenated Δ HbO and deoxygenated Δ HbR hemoglobin concentrations, denoising, and artifact removal. Namely, fNIRS data were firstly preprocessed to remove channels with a coefficient of variation greater than 10%; raw intensity signals were converted into optical density variations, corrected for motion artifacts using the temporal derivative distribution repair (TDDR) algorithm [35], and bandpass-filtered at (0.01-0.08) Hz.
The fNIRS protocol required 40 to 45 min on average, including the placement of the optodes, quality control of channels, and signal acquisition during the considered task, while the MRI protocol average 15 min required on. The pre-processing of fNIRS and MRI data were not time consuming. Conversely, a significant computational burden is referred to-the sensitivity matrix computation over SSA (Section 2.2). Nonetheless, the overall experimental procedure is suitable for clinical applications since the computation of the sensitivity matrix can be performed offline after a routine anatomical MRI is acquired. Moreover, the sensitivity matrix can be reused for several fNIRS trials in a follow-up protocol in most cases where a stationary anatomical condition can be assumed. The functional task consisted of a block-designed motor-grasping task, where participants had to alternatively move their left or right hand [8]. The fNIRS signals were preprocessed using the Nirstorm package (https://github.com/Nirstorm/nirstorm (accessed on 5 February 2023)) in Brainstorm [34] to extract the relative changes in oxygenated ∆[HbO 2 ] and deoxygenated ∆[HbR] hemoglobin concentrations, denoising, and artifact removal. Namely, fNIRS data were firstly preprocessed to remove channels with a coefficient of variation greater than 10%; raw intensity signals were converted into optical density variations, corrected for motion artifacts using the temporal derivative distribution repair (TDDR) algorithm [35], and bandpass-filtered at (0.01-0.08) Hz.
The fNIRS protocol required 40 to 45 min on average, including the placement of the optodes, quality control of channels, and signal acquisition during the considered task, while the MRI protocol average 15 min required on. The pre-processing of fNIRS and MRI data were not time consuming. Conversely, a significant computational burden is referred to-the sensitivity matrix computation over SSA (Section 2.2). Nonetheless, the overall experimental procedure is suitable for clinical applications since the computation of the sensitivity matrix can be performed offline after a routine anatomical MRI is acquired. Moreover, the sensitivity matrix can be reused for several fNIRS trials in a follow-up protocol in most cases where a stationary anatomical condition can be assumed.

Sensitivity Matrix Computation
Anatomical MRI data were processed in Nirstorm to compute the forward model at both wavelengths λ: Matrix A (λ) elements represent the cortical surface-to-channel sensitivity coefficients at both wavelengths. The channel vector of optical density data is ∆OD (λ) = [y i ], with i = 1, 2, . . . , I and I = 102 number of channels in our configuration. The vector α,j , j = 1, 2, . . . , J represents the absorption coefficient changes associated to the j th cortical element CE j activity, which in turn is related to hemoglobin concentration changes as ∆µ , where the α are the molar absorption coefficients of the respective chromophores at wavelength λ. The solution of this linear algebraic system is hence transferred from the surface channel level to the level of cortical elements CE j . Importantly, the CE j are automatically defined in Freesurfer by a vectorized mesh representation of about 15,000 elements for both subject-specific and atlas-based anatomies.
The optodes were placed on the scalp at the subject's MRI-native coordinates and a virtual 10/5 cap model was co-registered to the specific anatomy, assuring the correspondence of the relevant landmarks. The resulting grid of sources (red) and detectors (green) is shown in Figure 1.
For the sole purpose of a quality check, all scalp locations and optode landmarks were registered to the Montreal Neurological Institute (MNI) space, thus observing the low dispersion of the clusters of corresponding optodes through subjects and atlases (see Supplementary Material for details).
Matrix A λ was estimated using Monte Carlo (MC) simulations in the native MRI space of SSAs, or in the atlas coordinates for ABAs. The forward problem was approached with GPU-accelerated MC simulations [36] with 5 × 10 7 photons per optode. A 5-layer head model was considered: scalp, skull, CSF, GM, and WM. The optical parameters shown in Table 1 were set from the work of Tak et al. [37] and Eggbrecht et al. [27]. Next, the channelwise sensitivity profile was obtained from the respective source and detector fluencies [38], hence interpolated onto the pial surface (i.e., interface between GM and CSF, from now on simply referred to as cortical surface) using a Voronoi-based method [39] and smoothed with a 2 mm FWHM gaussian kernel. The same analysis was applied to three atlases. This ABA approach had the dual aim of assessing the validity of applications in which the absence of an individual SSA is surrogated by an ABA [19,33] and also to provide three benchmarks to our SSA group dispersion analysis. The selected atlases were Colin27 [40] due to its wide adoption in fNIRS applications [14] and software [7,41,42], ICBM152 [43] as a standard MNI template, and FSAverage as the default brain template implemented in the Freesurfer package [44].
A qualitative view was provided by mapping the total sensitivity (TS) of all channels vs. cortical elements CE j . This was readily computed by summing the values in each column of matrix A: The inverse problem (i.e., projecting surface measures onto the cortex) is ill-conditioned and undetermined, thus requiring heavy regularization. For the given application examples (see Section 3.4), we employed the solution based on the minimum norm estimate approach (for both λ, dropped for clarity), according to the work of Machado et al. [24]: where ĥ 1 ,ĥ 2 are regularization hyperparameters estimated through the restricted maximum likelihood method, while (C 1 , C 2 ) are the covariance matrices, respectively, employed to model measurement noise and the a priori distribution of absorption changes. Notably, Equation (1) indicates the size of regularization to solve for the ill-conditioned and undetermined problem presented in Equation (3), since estimating an order of 10.000 unknown cortical elements' values from an order of 100 surface measures. Conversely, the present study focuses the forward sensitivity matrix A itself, based on the sole SSA (or atlas) and the simulated optical properties. Namely, two aspects are considered: (i) the drop in sensitivity with lateral displacement, statistically analyzed across channels and subjects by sensitivity displacement surfaces (SDS, Section 2.3); (ii) the variability of sensitivity coefficients grouped by cortical areas into area to channel (A2Ch) sensitivities (Section 2.4).

Sensitivity Displacement Surfaces
This section describes the proposed method to statistically assess the mean 3D profile of sensitivity with respect to displacements from the channel center (CC) at increasing cortical depths. For this purpose, the longitudinal γ (i.e., in the source to detector direction) and transverse τ (i.e., orthogonal direction to γ) curvilinear distances from the CC were defined for each i th channel considering the scalp curvature specific to each subject. The depth d was defined as the distance from a given (γ, τ) scalp location along the local normal to the scalp down to a cortex element CE j , for which sensitivity A i,j was also recorded. Since the cortical element index j is univocally determined by the channel i and the displacement (γ, τ), the sampled A i,j can be written as A i (γ, τ). Hence, in a given subject (sbj) and in the frame of the i th -channel, a depth map d sbj i (γ, τ) and a sensitivity map A sbj i (γ, τ) were derived. Coordinates (γ, τ) were sampled over a 40 × 40 mm 2 field by 2 × 2 mm 2 bins (with some approximation at borders, due to the curvilinear head shape).
To gain a statistical representation of the sensitivity profile A sbj i (γ, τ, d), depth and sensitivity values from all 102 channels were superimposed by aligning the respective frames. Thus, for each (γ, τ) coordinate the relationship between depth and sensitivity could be extracted from 102 samples. A segmentation by depth was fixed according to the quartiles of the overall sampled depths in the subject: d sbj Q , with the quartile index Q = [1; 2; 3; 4] ranging from the least cortical depths (i.e., gyri) to the deepest ones (i.e., sulci). Thus, it was possible to map the average sensitivity profiles at progressive depths according to four sensitivity displacement Surfaces (SDS) averaged across all channels for a single subject: To quantify the expected drop in sensitivity with depth, the integral value of the four SDS, named integral under the surface (IUS), were computed and expressed in percentage values: where Q = [1; 2; 3; 4] indicates the quartile of cortical depth. The shape of the SDS curves was assessed by the full width at half-maximum (FWHM) along the longitudinal FWHM λ and transverse FWHM τ directions. Namely, FWHM λ and FWHM τ were computed by considering the profile of SDS curves along these directions and calculating the distance between points at half of maximum value. Values for the first and second quartiles of cortical depth are reported, since almost flat SDSs were found at further depths.

Area to Channel Sensitivity
The interpretation of surface fNIRS data are mainly addressed to detect the activations of cortical areas expected to be involved during the execution of a task or in response to specific stimuli. In this perspective, the cortical elements belonging to each area may be grouped as subset of cortical elements CE j composing an area a (i.e., cortical parcellation), which are identified by the subset of indexes J a ≡ j CE j ⊂ area a , with a = 1, 2 . . . N areas . Accordingly, the columns of A λ relevant to the same area a are summed up into the a th -column of a reduced sensitivity matrix B λ : The coefficients B λ i,a express the value of area to channel (A2Ch) sensitivities. This reduced model is adopted here to evaluate the sensitivity profile of a parcellation area a due to the specific i th -channel and across SSA vs. ABA cases.
The final A2Ch values were normalized using the sum of channel-wise sensitivity coefficients across cortical elements and expressed as percentage ratio.
As results, A2Ch%(i, a) values span in the [0, 1] range and express the relative influence of an area over a measurement channel (i.e., area to channel coupling) and its variability through subjects.
The A2Ch% values were analyzed on a representative set of Brodmann areas (BA) (Figure 2A-D panels): bilateral BA1-2-3-4 as primary sensorimotor cortex (SensoriMotor Network, SMN), BA17 as primary (Visual Network 1, VIS1), and BA18-19 as secondary visual areas of occipital cortex (Visual Network 2, VIS2), BA46-9 as high-level, and BA45-47 as low-level sites of executive functions of dorsolateral-(PFC1) and ventrolateralprefrontal cortex (PFC2), respectively. The correspondence between BAs with anatomical areas, along with the notation employed in the results section, is presented in Table 2.  We also considered the bilateral precentral and postcentral gyri as additional anatomical parcellations on the Destrieux atlas [45] ( Figure 2E-G panels). Separate analyses of the selected BAs and Destrieux areas were performed. In both cases, A2Ch% values above a threshold of 0.2 are presented and considered as significant channel-wise values (i.e., at least of 20% of the sensitivity profile of the considered area is explained by that channel), while empty channels (i.e., those distant from any considered areas) are omitted. Table 2. Summary of employed left (L-) and right (R-) hemisphere areas and their correspondence between Brodmann areas and anatomical areas.

Functional Label
Functional Areas Anatomical Areas Figure 3 provides an example of sensitivity profile of channel S9D7 (i.e., the representation of A i,j i=S9D7 at source 9 and detector 7), while the red dot represents the S9D7 midpoint as channel position. Considering the log 10 scale down to 1/100 (i.e., blue in the color scale), it is evident that the sensitivity profile drops in a short range of depth from the top of gyri to the bottom of sulci. Moreover, there is a noticeable change in the sensed pattern using SSAs (A-C panels) and ABAs (D-F panels), whose variations can be better-quantified with A2Ch% indexes (Section 3.3).  The total sensitivity maps on the cortical surface are shown in Figure 4 for three subjects and the three atlases. Values are normalized to the maximum sensitivity (i.e., red color scale) and in log • scale down to 100-fold less (i.e., blue color scale). The qualitative inspection of these maps anticipates the quantitative results of the following analyses. The most evident feature is that the drop in sensitivity with depth shown above for a single channel (S9D7, as example) is maintained, also summing up all channel sensitivities. This confirms that the signal recorded on the scalp surface is almost entirely determined by the gyri. Secondly, along the sagittal plane, the inflation shows a part of the The total sensitivity TS j maps on the cortical surface are shown in Figure 4 for three subjects and the three atlases. Values are normalized to the maximum sensitivity (i.e., red color scale) and in log 10 (·) scale down to 100-fold less (i.e., blue color scale). The qualitative inspection of these maps anticipates the quantitative results of the following analyses. The most evident feature is that the drop in sensitivity with depth shown above for a single channel (S9D7, as example) is maintained, also summing up all channel sensitivities. This confirms that the signal recorded on the scalp surface is almost entirely determined by the gyri. Secondly, along the sagittal plane, the inflation shows a part of the mesial cortex, which is a gray color scale since it is below the considered sensitivity range, hence negligible. Finally, similarities and differences among SSAs (A-C panels) and among atlases (D-panels) are highlighted. Overall, the maps presented in Figures 3 and 4 attain to MC simulation at λ = 760 nm; however, λ = 850 nm provided similar qualitative information.  Table 3 provides a summary of the percentage IUS values (median and IQR) at increasing quartiles of cortical depth. Results quantify how the sensitivity profile drops with depth. Additionally, Table 4 presents the depth separation values between each subsequent quartile: / , / , and / , respectively. The IUS reflects the integral of the sensitivity profile to cortical anatomy at the top of gyri and shows that the sensitivity profile is mostly confined to this cortical range. Values for SSAs had a median of IUS 0.39, while lower values are found for most ABAs: 0.3 for Colin27, 0.36 for ICBM152. Conversely, the IUS over the sensitivity profile of FSAverage is comparable to the one of SSAs.   Table 3 provides a summary of the percentage IUS values (median and IQR) at increasing quartiles of cortical depth. Results quantify how the sensitivity profile drops with depth. Additionally, Table 4 presents the depth separation values between each subsequent quartile: d Q1/Q2 , d Q2/Q3 , and d Q3/Q4 , respectively. The IUS 1 reflects the integral of the sensitivity profile to cortical anatomy at the top of gyri and shows that the sensitivity profile is mostly confined to this cortical range. Values for SSAs had a median of IUS 1 = 0.39, while lower values are found for most ABAs: 0.3 for Colin27, 0.36 for ICBM152. Conversely, the IUS 1 over the sensitivity profile of FSAverage is comparable to the one of SSAs.

Sensitivity Displacement Surfaces
Considering the sum of the two upper quartiles IUS 1 + IUS 2 values of 0.683 were found for SSAs, 0.596 for Colin27, 0.656 for ICBM152 and 0.681 for FSAaverage, which further confirms that the top of the sulcogyral pattern should be considered in the interpretation of fNIRS (i.e., almost of 70% of the sensitivity profile is entirely confined within this cortical range). Moreover, Table 4 highlights the steep drop of sensitivity with cortical depth. The quartile separation values in SSAs present a median decrease of about 2 mm per quartile (i.e., 11.8 mm at d Q1/Q2 , 13.6 mm at d Q2/Q3 , 15.7 mm at d Q3/Q4 ). Considering the shape of the SDS plots, it is possible to evaluate the longitudinal and transversal sensitivity decay. Figure 5 displays this analysis relevant to two SSAs. The drop in amplitude with the depth quartile (from SDS 1 to SDS 4 ), which is quantified by the IUS% indexes, is visually confirmed. The sensitivity profile with displacement reveals two major features: (i) high sensitivity is kept for a wide portion of the longitudinal γ displacement, yet with a drop approaching the source and the detector point; (ii) a sensibly steeper drop is seen in the transverse direction τ.  Considering the sum of the two upper quartiles IUS + IUS values of 0.683 we found for SSAs, 0.596 for Colin27, 0.656 for ICBM152 and 0.681 for FSAaverage, whi further confirms that the top of the sulcogyral pattern should be considered in the inte pretation of fNIRS (i.e., almost of 70% of the sensitivity profile is entirely confined with this cortical range). Moreover, Table 4 highlights the steep drop of sensitivity with cortic depth. The quartile separation values in SSAs present a median decrease of about 2 m per quartile (i.e., 11.8 mm at / , 13.6 mm at / , 15.7 mm at / ). Considering the shape of the SDS plots, it is possible to evaluate the longitudinal an transversal sensitivity decay. Figure 5 displays this analysis relevant to two SSAs. T drop in amplitude with the depth quartile (from SDS to SDS ), which is quantified the IUS% indexes, is visually confirmed. The sensitivity profile with displacement revea two major features: (i) high sensitivity is kept for a wide portion of the longitudinal displacement, yet with a drop approaching the source and the detector point; (ii) a sen bly steeper drop is seen in the transverse direction τ. These features were common to all SSAs and are quantified in Table 5   These features were common to all SSAs and are quantified in Table 5 by the full width at half-maximum (FWHM) along the longitudinal FWHM γ and transverse FWHM τ direc-tions, respectively. Values for SDS 1 and SDS 2 are reported since they justify most of the sensitivity. The remarks advanced about Figure 5 are confirmed: (i) FWHM γ covers about 21 mm of the longitudinal axis, which is about 60% of the source-detector average distance; (ii) the sensitivity profile of SSAs is sensibly confined along the longitudinal direction by the steeper transverse drop with FWHM τ about threefold narrower than FWHM y . Conversely, this result is not valid for the ABAs case, FWHM τ,1 and FWHM τ,2 being almost comparable or even higher than FWHM γ,1 and FWHM γ,2 , respectively. Therefore, this result suggests that ABAs provide a more uniform drop in the sensitivity profile along longitudinal and transverse direction, possibly due to the smoother anatomy resulting from merging many subjects.

Area to Channel Sensitivity
The A2Ch% coefficients are evaluated here, considering the sample areas defined in Section 2.4. Variability is presented by the SSA group statistics (median (IQR) and compared to the values of the three ABAs. Since the result shown in Section 3.2 displayed negligible differences in sensitivities at λ = 760 nm and 850 nm, the following analyses will be shown only for the former wavelength. Tables 6 and 7 report the results of A2Ch%coefficients for the selected Destrieux areas, while Tables 8-12 report functional ROIs associated with selected BAs. In both cases, results for SSAs are reported as median (IQR) values across subjects.
To limit the size of these tables, channels displaying an A2Ch% < 0.2 in all columns (i.e., SSA median and ABA values) were omitted. Despite this cutoff, the tables show that the considered areas are sensed by a fairly long list of channels, generally around 10. This indicates that the contribution of many channels enters with varying weights the quantification of the activation of a specific cortical area provided by the inverse problem.
Remarkably, in several areas, one or more channels were found with a median A2Ch% close to or above 0.5, hence indicating that those channels prevalently sensed that area. However, in many of these cases, a wide IQR was found. Hence, interpreting that channel as primarily indicative of that area is only worth it for some of the subjects.
SSA results were also compared to ABAs. The tables highlight A2Ch% values for which the SSA median differed by more than 0.2 from any one of the ABA cases in bold. Such differences, when largely present, indicate that surrogating the SSA with an ABA might provide misleading results when mapping surface data to cortical anatomy. Namely, the highest rate of SSA vs. ABA differences (i.e., highest number of highlighted table rows) is found around the VIS regions (Tables 8 and 9), while the lowest rate around the PFC regions (Tables 11 and 12). An intermediate case was represented by motor and sensorimotor regions, as highlighted in Tables 6 and 7 for Destrieux areas and Table 10 for the SMN region. Although the analysis was limited to few representative ROIs, it is worth noting that the IQR range of A2Ch% coefficients across SSAs is often comparable or higher than the median value, hence indicating that there is an intrinsic variability regarding the identification of cortical anatomy.     Figure 6 provides an example of cortical image reconstruction of the block-averaged precentral and postcentral gyri response to the right-hand motor task in a representative subject (i.e., subject #9). Surface measures, averaged from 10 repetitions, where projected according to the own SSA (A-C panels) and compared with those resulting from employing the three considered atlases: Colin27 (D-F panels), ICBM152 (G-I panels), and FSAaverage (J-L panels). The color maps represent the peak ∆[HbO 2 ] and ∆[HbR] values (left and mid column, respectively). The averaged responses are shown in the right column relevant to the contralateral precentral and postcentral gyri (integral in the areas). As expected, the time course shapes were negligibly affected by the adopted anatomy, since they are bound to the task paradigm and to the hemodynamic response function. Conversely, large differences were seen in their amplitudes, in keeping with the differences in the cortical maps. Namely, the peak of ∆[HbO 2 ] response (650 µmol/L) was found in the precentral gyrus if the correct SSA was applied. Conversely, the main activation was shifted on the postcentral gyrus for the ABAs, with large amplitude differences: 450 µmol/L for Colin27, 980 µmol/L for ICBM152, and 700 µmol/L for FSAverage. Importantly, the ratio between the postcentral and the precentral gyrus shows sensible differences, being about 2:1 Colin27 and FSAverage, while up to 3:1 for ICBM152.

Example of Cortical Image Reconstruction
This result confirms that the anatomy considered in the inverse problem solution strongly influences the localization and quantification of functional responses, which is explained by the forward problem analyses. Namely, the above results relevant to channel S9D7 can be revisited under this perspective. Qualitatively, in Figure 3, by the sensitivity patterns on different anatomies. Quantitively, by the variable A2Ch% values for the Destrieux left precentral and postcentral gyri in Table 6 and for the L-SMN in Table 10.
(left and mid column, respectively). The averaged responses are shown in the right column relevant to the contralateral precentral and postcentral gyri (integral in the areas). As expected, the time course shapes were negligibly affected by the adopted anatomy, since they are bound to the task paradigm and to the hemodynamic response function. Conversely, large differences were seen in their amplitudes, in keeping with the differences in the cortical maps. Namely, the peak of Δ[HbO ] response (650 μmol/L) was found in the precentral gyrus if the correct SSA was applied. Conversely, the main activation was shifted on the postcentral gyrus for the ABAs, with large amplitude differences: 450 μmol/L for Colin27, 980 μmol/L for ICBM152, and 700 μmol/L for FSAverage. Importantly, the ratio between the postcentral and the precentral gyrus shows sensible differences, being about 2:1 Colin27 and FSAverage, while up to 3:1 for ICBM152.
This result confirms that the anatomy considered in the inverse problem solution strongly influences the localization and quantification of functional responses, which is explained by the forward problem analyses. Namely, the above results relevant to channel S9D7 can be revisited under this perspective. Qualitatively, in Figure 3, by the sensitivity patterns on different anatomies. Quantitively, by the variable A2Ch% values for the Destrieux left precentral and postcentral gyri in Table 6 and for the L-SMN in Table 10.

Discussion
In this work, we showed an approach to (i) statistically assess the distribution of fNIRS sensitivity profile with respect to channel displacements and scalp to cortex distance; (ii) quantify the coupling between channels and selected functional/anatomical cortical areas (i.e., A2Ch% sensitivity). The proposed method addressed the multimodal integration of fNIRS and MRI techniques to further support the applicability of fNIRS in clinical research by alleviating the problems related to its limited spatial resolution and indirect sensing of cortical areas. Indeed, clinical applications would require a specific identification of cortical activations in patients with impaired neurovascular activity. The pairing of fNIRS with other recognized clinical neuroimaging techniques, such as anatomical MRI, is essential for progress in this direction while simultaneously taking advantage of its portability and potential applicability in various clinical contexts. For these reasons, our analysis was limited to multichannel fNIRS equipment with optodes placed at 10/5 locations, while the ongoing progresses relevant to more sophisticated instrumentation, optimized optode placement, and/or high-density probes were out of our scope.
Results for SDSs and IUS values provided a numerical and visual outcome regarding the spatial sensitivity distribution of the employed fNIRS probes along longitudinal and transverse direction at varying cortical depths. Results showed the characteristic features of the sensing field in SSAs: (i) along the longitudinal axis (i.e., source-detector direction) the sensitivity profile preserved meaningful values along a high percentage of the entire source-detector distance (i.e., the FWHM γ is about 21 mm across both the first and second quartiles of depth), while along the transverse axis we found a steeper decay in sensitivity (i.e., FWHM τ is about 6 mm across both the first and second quartiles of depth); (ii) a very wide drop in sensitivity is seen with depth, such that the bottom of sulci was attenuated about 100fold compared to the top of the gyri. This last result is quantified by the integral values of SDS plots according to quartiles of cortical depth. Approximately, the first depth quartile (i.e., top of gyri) provided about 40% of the sensitivity, the second about 30%, the third about 20%, and the fourth (i.e., bottom of sulci) about 10%. This result showed limited group dispersion across the SSAs, as quantified by low IQRs, since the SDSs were obtained by combining all channels, thus compensating the individual differences in the patterns of gyri and sulci. However, it should be stressed that this result refers to a group of healthy young adults with similar age span. Individual variability of fNIRS sensing with depth may conversely show even wider dispersions incase significant age spans were considered, or even neuropathologies implying brain atrophy [46][47][48].
The comparison of IUS values with the ABAs cases delivered similar results except for the Colin27 atlas, which provided an underestimation of the first quartile. In general, the proposed analysis is in line with recent works. Namely, Tian and Liu [49] proposed ICBM152 as reference for a depth compensation algorithm to improve the localization of functional activation. Liu et al. [50] adopted SSAs, indicating that fNIRS surface measurements can infer deep-brain activity comparable to simultaneous fMRI acquisitions. Strangman et al. [14,51] mapped the depth sensitivity profile over a Colin27 function of source-detector separation using MC simulations and quantified the influence of scalp and skull thickness. Moreover, the interpretation of fNIRS findings can be highly affected by anatomical variability and optode placement [15].
Regarding the A2Ch% analysis, in the fNIRS literature, it is a common practice to constrain scalp measurements to the cortical surface [21,33]. However, to the best of our knowledge, a numerical assessment of the forward problem and individual anatomic variability is usually overlooked, while directly assessing the inverse problem solution and regularization. The analysis was necessarily limited to a few Brodmann areas (BAs), yet representative of frontal, sensorimotor, and occipital regions. Separately, the precentral and postcentral gyri from the Destrieux atlas were also considered, given the different logic of this atlas based on anatomical features. To summarize the huge and sparse intersection between areas and channels, we provided a list of channels for each area, showing A2Ch% ≥ 20%. The major result was the complex portrait of channels sensing each of the considered areas that provided information about the area activation while solving the inverse problem. Importantly, the SSA group analysis revealed high IQRs through the subjects. Finally, for many A2Ch%, the result provided by one or more of the three atlases sensibly differed from the SSA group median. The major deviation of the ABAs from the SSA median were found in the occipital visual regions, decreased in the sensorimotor region, and even more in the frontal region.
The results of this study show that the further computational effort needed for cortical fNIRS mapping may correct the nonnegligible individual variability of the sensitivity of channels vs. cortical areas, thus improving interpretation via the inclusion of subjectspecific anatomy. Indeed, future developments of this work will addressed the inverse problem of providing a more stable solution by reducing the number of sensitivity matrix coefficients. Namely, the proposed method will allow for the definition of a spatial prior on cortical depth to avoid the inversion of sensitivity matrix elements not significantly contributing to the sensitivity profile (i.e., the SDS analysis of Section 2.3) and/or that do not properly map a cortical area of interest (i.e., the A2Ch% analysis of Section 2.4).

Conclusions
We reconsidered the current problem of fNIRS image reconstruction by introducing methods for the numerical quantification of the sensitivity profile with respect to cortical depth and coupling with cortical areas of interest. This step requires the estimation of a sensitivity matrix of surface channels vs. cortical elements, whose dimension challenges qualitative and quantitative analyses. We also limited the analysis to a multichannel fNIRS equipment according to a 10/5 system configuration of optodes. The choice was derived from the need for providing concrete indications for the clinical research use of fNIRS. An indication was obtained using an analysis of the cortical surface sensed on average by each scalp channel, which is limited in the transverse direction (orthogonal to the sourcedetector line) and drops dramatically with depth. The sharp sensing volume and the subjective variability of gyri and sulci, even in a group of young healthy adults, provided high dispersion of the coupling between channels and cortical areas of major interest, as indicated by the area of channel sensitivity coefficients. The analysis performed on three major atlases showed that such standard anatomies may provide a limited surrogate for subject-specific anatomies according to different regions of interest. In conclusion, moving towards the application of fNIRS in clinical contexts would greatly benefit from the integration of fNIRS data with the subject-specific MRI anatomy as routine practice.