Continuous wave functional near-infrared spectroscopy (CW-fNIRS) is a noninvasive imaging technique capable of measuring and mapping local changes in cerebral oxygenation caused by brain metabolic and hemodynamic activity.1 In CW-fNIRS, a set of optic fibers, called optodes, is placed on the surface of the scalp, emitting (sources) and collecting (detectors) NIR light (16 dual sources and 32 detectors with our Brainsight CW system). Biological tissues are highly scattering but present low absorption in the NIR optical window (650 to 950 nm). Therefore, a portion of emitted photons can travel through the head up to a depth of 2 to 3 cm within the tissues following a diffusive process, reach the outermost cortex, and then be backscattered to the scalp and eventually to the detectors.2 In fNIRS, the measurement channels on the scalp are defined by the different source/detector pairs (SD pairs) within a certain separation range. Indeed, the larger the separation between a source and a detector, the greater the average depth of penetration of the photons in the brain and the lower the measured diffuse reflectance signal. Based on theoretical investigations, it was demonstrated that a mean SD separation in the range of 2 to 4 cm is a good trade-off between penetration depth and signal-to-noise ratio (SNR).3
Brain activity leads to an increase in oxygen consumption, which is accompanied by an increase in cerebral blood flow due to neurovascular coupling.4 This brings about changes in local oxygenated (HbO) and deoxygenated (HbR) hemoglobin concentrations and, therefore, local variations in the tissue absorption properties. HbO and HbR are two dominant NIR chromophores in cerebral tissues. Consequently, during neuronal activation periods, the variations in the amount of diffuse reflectance measured over a set of SD pairs on the scalp can be used to compute the associated changes in [HbO] and [HbR] with respect to a baseline (in what follows, [X] will denote the concentration of X and will denote a variation relative to a baseline). Measurements at multiple wavelengths (e.g., 690 and 830 nm) are required to extract the concentration of each chromophore.
Either the tomographic imaging approach denoted by diffuse optical tomography (DOT)5 or the topographic imaging approach denoted by functional near-infrared topography (fNIRT)6,7 may be used in fNIRS to convert and map the variations of optical densities () measured on the scalp into and activity. DOT reconstructs and activity in the volume by solving an ill-posed inverse problem,8,9 whereas fNIRT produces topographic maps on the scalp surface using the differential modified Beer–Lambert law.10,11
In fNIRS, a specific arrangement of sources and detectors on the scalp is called an SD montage and most fNIRS acquisitions are conducted using semirigid pads with the optodes arranged in a grid mesh. Early fNIRS studies were conducted with sparse, single-distance imaging arrays with the difference in locations of adjacent sources and detectors being , in order to obtain measurements of similar and significant cortical sensitivity. Topographic and tomographic images obtained with these sparse arrangements, however, were usually mislocalized, blurred, and distorted with respect to their original shape because of the diffuse nature of light propagation and the sparse sampling provided by the SD montages.12,13 It has been demonstrated using phantoms, numerical simulations, and in vivo experiments that the spatial resolution of images can be considerably improved using dense SD montages with overlapping measurements at multiple SD pair separations. These results have been established for fNIRT12,14184.108.40.206.–19 and for DOT.13,2021.22.23.–24
Because it generally creates no interference with other imaging devices, fNIRS can be used in conjunction with other modalities, such as electroencephalography (EEG),25 functional magnetic resonance imaging (fMRI),26,27 and magnetoencephalography (MEG).28 Therefore, fNIRS has become an attractive option for the development of clinical applications that investigate neurocognitive processes associated with brain disorders.29 In epilepsy research, simultaneous EEG/fNIRS offers an ideal multimodal approach for studying the relation between hemodynamic activity and spontaneous epileptic discharges detected on the scalp EEG.3031.32.33.–34 It is important to find an optimal arrangement of both the electrodes and the available sources and detectors on the scalp in order to obtain an accurate estimation of the underlying hemodynamic activity; this issue is especially important because the presumed epileptogenic focus under investigation is patient-specific and may involve different brain regions. The combination of EEG and fNIRS in the context of epilepsy, however, entails constraints on the design of SD montages.
First, in EEG/fNIRS, one does not generally use an array of optodes embedded into semirigid pads because the placing constraints on the patches conflict with the EEG montages. The pads have to be strapped on the head around the the targeted region and the arrangement of the optodes has constant separations; on the other hand, the EEG electrodes must be positioned using a relative positioning system, such as the international 10/20 system or 10/10 system.35,36 EEG/fNIRS caps or rigid whole-head helmets allow greater flexibility in positioning optodes.37 They allow EEG recordings and provide a set of possible positions (called optode holders) for holding the source or detector optodes. Contrary to what is the case in semirigid pads, the separations between the nearest and next nearest optode holders on EEG/fNIRS caps may span a large range of values and may vary by several millimeters from one subject to another because EEG/fNIRS caps are generally made of a stretchable material in order to accommodate the scalp curvature. Using EEG locations for optode holders allows for standardization of optode placement across subjects and should ensure that the brain region explored by a particular SD pair is the same across subjects taking into account the variability specific of the anatomy.38 Even so, since the sensitivity of fNIRS measurements depends strongly on the thickness of the superficial tissue,39 it is difficult to ensure, when using either constant or variable SD separation measurements, that a given cortical area can be measured by several SD pairs with an appropriate and homogeneous SNR level.12,39 With prior information on optical sensitivity, however, the differences in holder separations on an EEG/fNIRS cap may be exploited to design adequate SD montages that are subject-specific. This is the strategy we considered in this study.
Second, many studies involving multiple subjects will use the same SD montage positioned approximately at the same location in order to perform group analysis. In EEG/fNIRS epilepsy studies, however, the epileptogenic regions to be explored with fNIRS are patient-specific.40 Typically, the patient-specific volumes of interest (VOIs) to be explored in priority with fNIRS may result, for instance, from the segmentation of an underlying lesion in anatomical MRI,41 from the source localization of EEG42,43 or MEG discharges,44,45 and from the fMRI blood oxygen level dependent response to epileptic discharges.46,47 Therefore, an SD arrangement has to be designed and optimized for each individual patient. Moreover, epilepsy is organized as an underlying epileptogenic network,48 and the simultaneous exploration of different VOIs may yield useful information. In that case, a limited number of sources and detectors must be distributed over the different VOIs.
In the present article, we describe an algorithm for computing the optimal set of positions on an EEG/fNIRS cap, where optimal means that the SD montage provides the best sensitivity to a target VOI. To do so, patient-specific light sensitivity profiles were estimated using Monte Carlo (MC) simulations and the optimization problem was formulated as a mixed linear integer programming problem.49
Materials and Methods
Isometric Cap and the 10/05 EEG/fNIRS Cap
To study the hemodynamic response associated with epileptic activity, an EEG/fNIRS cap is installed on the head of the patient. The cap must fit various head shapes and offers great flexibility for SD arrangements, which is of paramount importance when designing an SD montage for a specific epilepsy patient. In this study, we propose to evaluate two possible electrode/optode holder arrangements.
The first is an isometric arrangement [Fig. 1(a)] that we have designed in collaboration with EasyCap (Germany) in order to perform EEG/fNIRS acquisitions with epilepsy patients. The cap can accommodate 63 EEG electrodes and up to 123 optodes all over the head. On this cap, both EEG electrodes [in black in Fig. 1(a)] and optode holders [in green in Fig. 1(a)] are distributed following an isometric arrangement in concentric circles. There is, on average, a distance of 3 cm between two concentric circles of optode holders; along one circle, the average distance between two holders is set to 1.5 cm. EEG circles are interleaved with the optode holder circles and do not follow a standard EEG positioning system.
Jurcak et al. have proposed to standardize optode locations on the scalp using derivatives of the international 10/20 EEG positioning system,38 such as the high-density 10/05 system.50 Therefore, the second arrangement was a 10/05 montage. EEG/fNIRS caps using a subset of this arrangement have already been proposed37 [see Fig. 1(c)]. In the present article, for validation purposes, we designed a complete virtual 10/05 EEG/fNIRS cap on the skin mesh of one subject [see Fig. 1(d) and Sec. 3], following the approach proposed in Perdue et al.39 The virtual 10/05 cap can accommodate 69 EEG electrodes (not shown) all over the head at the EEG 10/10 system positions and up to 248 optode holders [the black points in Fig. 1(d)] at the remaining positions of the EEG 10/05 system.
NIR Light Propagation in Biological Tissues: Monte Carlo Simulations
Given that in human tissue scattering generally dominates absorption, and assuming that any variation in light intensity occurs on a time-scale much greater than the mean transit time of light, the light propagation can be modeled by the time-independent diffusion equation.51
is the photon fluence rate () at wavelength . is the position vector in the volume space . is the photon diffusion coefficient with the reduced scattering coefficient (), which is the reciprocal of the photon random walk step length. is the absorption coefficient (). is the source distribution of photons ().
During brain activity, regional changes in blood flow and volume alter the concentration of hemoglobin in the brain and, therefore, change the absorption of the tissue.
denotes the baseline absorption. If we consider an optical measurement between a source and detector positioned at and , the solution of the photon diffusion equation after the perturbation is obtained using the Rytov approximation.51
is the fluence distribution at the baseline state (before the perturbation). is the Green’s function of the diffusion equation at the baseline state. The Green’s function is the solution of Eq. (1) considering an infinitely small point source (represented analytically as a Dirac delta function in space) at position . The ratio of light that reaches the detector after a variation in the underlying tissue absorption is reported as a change in optical density () and defined as
Using Eq. (5) in a discrete volume ( voxels), the optical density can be rewritten as follows:6) can be rewritten using the following matrix product:
is a column vector describing the change in OD for each measurement , is a column vector describing the change in absorption for each voxel in the volume, and is a matrix (number of measurements by number of voxels) whose each row is a vector containing the scalar sensitivity coefficients (mm) of a particular SD measurement to changes in absorptions within each voxel . In , the solutions and can be computed using analytical or numerical methods of light propagation in biological tissues.8 For realistic optical montages and complex head geometries, the numerical solutions methods, such as finite element or MC simulation, are mandatory.51,52 Indeed, the sensitivity profile depends on the local individual anatomy.3,53,54
In this study, the solution associated with each optode holder has been obtained through an MC approach,52 which can take into account the specific anatomy of each patient using a realistic multilayered head model built from the anatomical magnetic resonance imaging (MRI) of the subject. The MC method simulates photon transport in tissues. A photon packet with a specific weight is launched into the medium from a source point consisting of the accurate position of one of the optode holders coregistered to the anatomical MRI. The interaction of the photon packet with matter is treated as a stochastic process, where each randomly selected interaction causes a random change in the direction of the path of the photon packet propagation as well as in its weight, until it is either terminated or reflected. Any number of photon packets can be launched and modeled, until the resulting solution has the desired statistical SNR. The MC method is equivalent to the modeling of photon transport analytically by the radiative transfer equation, and it provides accurate solutions especially in the cerebrospinal fluid (CSF) (where the diffusion approximation may not hold). To decrease the computational burden, we used the MC method implemented in the Monte Carlo eXtreme software by Fang and Boas.55 This implementation uses massively parallel computing techniques, thanks to graphics processing units (2 Nvidia Tesla cards), and is extremely fast () when compared with traditional single-threaded CPU-based simulations.
Optimization of the SD Montage over the VOIs
The goal of the present study is to propose a method for finding the discrete positions of sources and detectors among the set of possible optode holder positions on an EEG/fNIRS cap. The choice of the discrete positions should maximize the spatial sensitivity profile of OD measurements in the target VOIs.
We formulated the problem as a mixed integer linear programming problem. Let (resp. ) denote the number of sources (resp. detectors) to be installed and the set of holder position indices. The variable (for ) is a binary variable that equals 1 if and only if a source is placed in the holder . Similarly, (for ) is a binary variable that equals 1 if and only if a detector is placed in the holder . For an SD measurement between holders and , denotes the sum of the sensitivity coefficients over the voxels () in given VOIs.6). In our study, we generated an optimal montage for both wavelengths, and was defined by using for each voxel the average of the and coefficients, which have similar orders of magnitude. Finally, we define
1. Any source or detector is positioned at exactly one optode holder, that is,
2. The definition of can be implemented by the following constraints:
(3) A source and a detector cannot be positioned at the same holder, that is,
This mixed integer linear programming problem was solved using a branch-and-bound algorithm; we refer the reader to Land and Doig49 and to the book by Nemhauser and Wolsey.56 In this case, the branch-and-bound algorithm consists of solving the linear relaxation of the above model, i.e., replacing the constraints that and are binary variables by the constraints that and are nonnegative real variables. Then the problem can be solved by the simplex algorithm for linear programming (see, for instance, Ref. 57). If the solution of the linear program includes a fractional variable (for instance, if equals 0.5 for some ), one creates two subproblems: one in which equals 0 and one in which equals 1. If both of these subproblems have integer solutions, they are compared and the best one is an optimal solution of the original problem. Otherwise new subproblems are created. The reader will find in the references quoted above a description of the whole procedure, including techniques for accelerating the branch-and-bound procedure. Algorithms for linear and integer linear programming have been implemented in several software. We have solved our model using MATLAB® (MathWorks, Massachusetts), which includes the ILOG CPLEX MATLAB® toolbox (IBM, version 13, http://www-01.ibm.com/software/commerce/optimization/cplex-optimizer/).
Isometric and 10/05 Arrangements on the Anatomical Head Model
In order to initialize the optical properties required for MC simulations, an anatomical head model was built from the classification into five tissue types [scalp, skull, CSF, gray matter (GM), and white matter (WM)] of the T1-MRI-scan of a human subject, previously acquired under an fMRI/fNIRS protocol approved by the ethics committee of the Montreal Neurological Institute and the McConnell Brain Imaging Center. The T1-MRI was obtained on a Siemens MAGNETOM 3T MRI scanner [repetition time , echo time ), ]. The MRI volume was contained within an array of with a resolution volume of . The classification approach5859.–60 combines an iterative, hierarchical, nonlinear registration procedure and a three-dimensional (3-D) digital model of human brain anatomy containing both volumetric intensity-based data and a geometric atlas.61 A slice from the classified head volume is shown in Fig. 2(a).
Based on estimates available from the literature,62 absorption and scattering coefficient values at 830 and 690 nm were associated with each tissue type constituting the segmented anatomical head model (Table 1). The anisotropy coefficient of tissues was set to 0.9 and the refractive indices of air and tissue were assumed to be 1.0 and 1.37, respectively.
Absorption/scattering (mm−1) coefficients used for the Monte Carlo simulations.
|Tissues||690 nm||830 nm|
After an appropriate placement of the isometric EEG/fNIRS cap on the subject head, the 123 holder positions were digitized using the Brainsight 3-D localizer device and coregistered to the anatomical MRI of the subject using a rigid body transformation63 based on anatomical landmarks (tip of the nose, nasion, left and right preauricular points). A mesh of the subject skin (86,111 vertices, mean distance between vertices: 1.4 mm ) was also obtained from the T1-MRI using the Brainstorm software.64 The 123 holder positions were projected on the skin mesh [Fig. 1(b)]. Considering only SD pairs with interoptode separations between 1 and 5 cm, this isometric optode holder arrangement allowed for 732 possible SD measurements.
In order to test the performance of the 10/05 arrangement, we also designed a virtual 10/05 arrangement with 317 positions on the same skin mesh, following the approach proposed by Perdue et al.39 Four landmarks (nasion, inion, and left and right preauricular points) were manually identified on the skin mesh. Then the NFRI MATLAB® toolbox38 was used to find the 10/10 positions on the head model. A linear transformation matrix was calculated from the 10/10 coordinate locations in Montréal Neurological Institute (MNI) space50 to the 10/10 coordinates on the scalp of the segmented head model. This transform was applied to the MNI 10/05 coordinate locations50 to find the remaining 10/05 locations on the model. For this study, the 10/10 positions (69 points) were reserved for EEG electrodes, whereas the remaining 248 positions were defined as potential candidates for fNIRS optodes [Fig. 1(d)]. Keeping only pairs with interoptode separations between 1 and 5 cm, this configuration allowed for a total of 2645 possible SD measurements.
Computation of Sensitivity Coefficients
Green’s functions of each optode holder (located at voxel ) of the isometric and 10/05 arrangements were estimated by MC simulations. For each simulation, 100 million photons with a unit survival weight were launched in order to provide results with low statistical noise. Eight (resp. 16) hours were necessary to compute the Green’s functions at 830 and 690 nm for all the holder positions of the isometric (resp. the 10/05) arrangement. The sensitivity profiles of each of the 732 (resp. 2645) possible pairs of optode holders on the isometric (resp. 10/05) arrangement were calculated using Eq. (6). A slice of one sensitivity profile at 830 nm for a particular SD measurement is shown in Fig. 2(b) and overlaid on the anatomical MRI of the subject.
Simulation of Regional Brain Activity in the Cortex
To verify that the proposed optimization method actually finds an adequate montage, we simulated 20 configurations with one focal (20 mm diameter) spherical VOI at different locations in the left hemisphere of the anatomical head model. In epilepsy research, it is often necessary to compare brain activity from the presumed focus with its corresponding homologous contralateral brain area, assumed to be normal. This is the reason why we simulated VOIs in the left hemisphere only. In such a configuration, fNIRS data acquisition will consist of mirroring the montage obtained on one side. The VOI centroids were positioned randomly in the GM with a minimal separation constraint of 10 mm between them and a maximal distance from the skin layer of 25 mm. The volume of each VOI was constrained within the GM and each voxel lying outside the GM was removed. We also simulated 20 configurations with one extended spherical VOI (50 mm diameter). The same centroids were used. Figures 3(a) and 3(b) illustrate examples of one focal versus one extended VOI.
To evaluate how the optimal positioning algorithm behaves when multiple VOIs need to be targeted simultaneously, we also simulated 20 configurations composed of two focal VOIs (with a minimum separation of 5 cm). Finally, we simulated 20 configurations involving two mirrored VOIs (i.e., a total of four VOIs). These configurations were proposed to mimic situations in which we plan to use fNIRS to target two different brain regions with their corresponding contralateral regions. We will refer to these last two configurations as the double VOIs and symmetrical VOIs configurations. Figures 3(c) and 3(d) illustrate, respectively, a double and a symmetrical configuration.
In each simulated VOI, we assumed homogeneous changes in absorption of and , mimicking [HbO] and [HbR] changes of from their baseline concentration; this is a realistic estimation of the variation usually encountered during functional activations. For any montage, the change in optical density at 830 and 690 nm can be calculated using Eq. (7) for each SD pair.
Parameters of the Optimization Model
Among all the 16 dual-wavelength sources and 32 detectors available with the Brainsight fNIRS system, we used half the optodes available, which would allow us to mirror the final montage to acquire homologous contralateral regions as well. Therefore, the configuration parameters of the optimization model are sources and detectors among a choice of optode holders when performing the optimization on the EEG/fNIRS cap with the isometric arrangement [Fig. 1(b)] and optode holders when performing the optimization on the virtual 10/05 arrangement [Fig. 1(d)]. In the simulation configuration involving two symmetrical VOIs, we used all the optodes available on the system, i.e., in this case, sources and detectors. The branch-and-bound algorithm was launched for each of the simulated VOI configurations and each optode holder arrangement separately.
Comparison of the Optimal Montages with Regular Arrangements of Sources and Detectors
Once the isometric optimal () and 10/05 optimal () montages were obtained for each simulated VOI, and were calculated for each SD pair of the montage. Then the topography was estimated by mapping the resulting value at each channel location (or sampling point). A channel location was defined as the middle point between the corresponding source and detector. The topography was then interpolated on the skin mesh using an inverse distance weighting interpolation scheme as in Aihara et al.65 and Takeuchi et al.66 The topographic maps were finally compared with the maps obtained when using three other arrangements of sources and detectors that we will consider as references.
1. Double density (DD) montages are often considered to provide an optimal contrast-to-noise ratio in activation imaging.12,14220.127.116.11.–19 Therefore, we compared our and montages with a virtual DD montage composed of 8 sources and 11 detectors (24 pairs). The optodes were interleaved in a specific pattern as shown in Fig. 4(a). For each VOI, the DD montage was positioned automatically over the volume to illuminate. To do so, the centroid of the simulated VOI was first orthogonally projected on the skin surface. We then designed a two-dimensional (2-D) planar DD montage centered on this point and oriented tangentially to the skin surface. The optode positions of this 2-D montage were finally projected on the skin surface. For all VOIs, the resulting DD grids were regular with a mean distance of 15 mm () between neighboring optodes. For all DD montages, MC simulations were performed to calculate the corresponding sensitivity matrices. Because a DD montage has a limited spatial coverage of the head, no DD montages were created for the double and symmetrical configurations.
2. With the isometric EEG/fNIRS cap, a simple configuration is to arrange sources and detectors in alternating circles at optode holder locations. Therefore, we designed a virtual whole-head regular isometric () montage on the EEG/fNIRS cap with 60 sources and 63 detectors (466 pairs) as illustrated in Fig. 4(b). Because we had already calculated the sensitivity profiles of each possible pair on the isometric EEG/fNIRS cap, no additional MC simulations were required. Obviously, such a whole-head configuration is not feasible experimentally with the typical number of sources and detectors currently available with fNIRS systems.
3. A virtual whole-head configuration following the 10/05 system with sources and detectors arranged in alternating coronal rows has been proposed by Perdue et al.,39 who demonstrated that this arrangement provides a relatively uniform sensitivity across the cortex. Using the available 248 optode positions of the virtual 10/05 cap, we designed a whole-head regular 10/05 () montage with 92 sources and 156 detectors (1250 pairs) arranged in alternating coronal rows [Fig. 4(c)]. Because we had already calculated the sensitivity profiles of every possible pair on the virtual 10/05 EEG/fNIRS cap, no additional MC simulations were required. Again, such a whole-head configuration is not feasible experimentally but will be used for comparison.
For each simulation configuration, a binary ground truth (GT) map on the skin was first defined as follows: we first identified the vertex of the skin surface with minimal distance to the center of mass of the VOI. All other voxels of the VOI were then projected to the skin surface following a normal to the tangential plane at this location. Let denote the number of active vertices () in this map denoted by . In the same way, let , , , , and denote the topographies of the optimal isometric, regular isometric, optimal 10/05, regular 10/05, and DD montages. The skin vertices set was then divided in two regions. As shown in Figs. 5(a) and 5(b), the first region illustrated in red was the projected VOI that we considered as the GT. The second region illustrated in green, called region of background, was defined as all connected vertices along the skin surface in the neighborhood of the GT up to an order of 15 connections along the mesh. will denote the number of vertices in the region of background.
To quantify the performances of the optimal montages when compared to the DD and whole-head regular montages for each VOI configuration, we propose three validation metrics. Validation metrics were evaluated using topographic maps obtained at 690 and 830 nm.
1. Detection accuracy: We propose a criterion based on receiver operating characteristic (ROC) curves67 following the approach suggested by Grova et al.68 Let denote the topography of one of the five montages under investigation at one wavelength. For a threshold chosen in the interval [0,1], we considered the ’th vertices of as active if its normalized energy is . Using the as a reference, we were able to quantify the sensitivity and specificity for each threshold and then draw the ROC curves. The area under the ROC curve (AUC) was used as detection accuracy index. It quantifies the detection ability in a close neighborhood of the GT, assessing whether the activity and its spatial extent were accurately recovered or not. Any overestimation or underestimation of the spatial extent of the GT should then affect this criterion. To interpret the AUC as a detection accuracy index, one should provide the same number of active and inactive vertices to the ROC analysis. Therefore, to calculate the amount of false positives and true negatives, we drew randomly vertices within the region of background. We performed this estimation randomly 20 times in order to average the AUC scores and get robust statistics (for further details see Grova et al.68).
2. Sensitivity: To assess the ability of the optimal montage to map the volumetric absorption changes into OD signals with maximum amplitude, an absolute local energy criterion (LE) was computed as follows:
3. Optimal sensing: To assess the amount of measurements overlap in each voxel, an overlap map was generated by thresholding and binarizing each row of the sensitivity matrix at 2% of its maximum value (as suggested in Tian et al.13). The number of overlapping measurements at each voxel was counted as the sum of contributions from all SD pairs for a given montage. The 2% threshold was chosen empirically after some tests have been performed in order to provide measurement density maps that have nonzeros coefficients within the cortical aspects of the GM volume. From this overlap map, we calculated the ratio of visible voxels (overlap score ) within the VOI and the mean value of the overlap score among the visible voxels.
The process of optimization was fast; for one montage the convergence of the branch-and-bound algorithm was reached in only few seconds on a processor IntelCore2 Quad CPU Q940.
For each montage type and VOI configuration (single focal, single extended, double, and symmetrical), Figs. 6 and 7 show the topographic maps generated by an absorption change in one simulated VOI. The corresponding sources (blue) and detectors (green) for each montage are shown on the scalp surface. The maps, represented using a yellow-to-red linear color scale, were normalized by their maximum value and the black isolines in each map indicate the half-maximum values. In this section, we present only results obtained at 830 nm since results at 690 nm were very similar.
For the simple focal VOI (Fig. 6, left), the maps of all montages are well localized with the maximum value lying in the GT area [Fig. 6(a)]. A statistical analysis performed on the 20 VOIs revealed that the mean Euclidean distance between the center of gravity of the GT maps and the estimated maps was for all montages. The topographic map associated with the DD montage [Fig. 6(b)] was the closest to the GT map. This could be expected since the DD montage was optimized to be exactly centered on the target VOI, whereas all other montages were optimized on the fixed optode holder positions. The and maps [Figs. 6(c) and 6(d)] were more extended than the GT area and slightly distorted. In this particular case, it is interesting to notice that the optimal positioning algorithm found positions for the eight sources, at exact positions of the regular isometric montage. This result has been observed for 3 out of 20 simulated VOIs. The and maps [Figs. 6(e) and 6(f)] appeared less blurred than the isometric montages likely due to the denser spatial sampling. The optimal positioning algorithm attributed five out of eight sources near the GT area. Three sources were positioned in the periphery and only provide weak amplitude.
For the extended VOI (Fig. 6, right), the DD map was again the closest to the GT map. and maps provided more spatially extended topographies than the GT. The optimal positioning algorithm found almost the same source positions for the focal versus the extended target VOI except for one source, which was moved to the middle row in order to capture extended activity in this region. The topography, which was irregular, was exhibiting a typical interpolation artifact of the inverse distance weighting scheme, caused by irregular spatial sampling. Indeed, the irregular distances between the optode holders in both the isometric and 10/05 arrangements lead to a nonhomogeneous density of measurement points over the scalp. This artifact did not appear with a DD montage where the sampling points were equidistant to each other. Whereas the was able to capture the full extent of the GT area, the maps missed some activity in the center of the GT map, probably because the corresponding SD measurements in the center had a low sensitivity to the targeted VOI.
When adding more than one spatial constraint to the optimization (example in Fig. 7), the double focal VOI and the double symmetrical VOI configurations, and montages were similarly sensitive to all the target VOIs. On the other hand, the complete montage was not able to recover activity from the most posterior VOI likely because the VOI was located too deeply in the brain for its sensitivity profile. Such behavior was observed for 2 out of the 20 simulated double VOIs with the montage.
For the (respectively the ), the branch-and-bound algorithm did not allocate positions for optodes over one the VOI for 1 out of 20 simulated double VOIs (respectively 5 out of 20). Such behavior of the algorithm was expected when working with an EEG/fNIRS cap, which features a high density of optode holders, such as the 10/05 cap. In such case, the optimization algorithm may favor one VOI (more likely a superficial one). A weight could be a priori attributed to each different VOI in the model to prevent this type of behavior.
Figure 8 displays the distribution of the AUC scores obtained over the 20 simulations and grouped by montage types. For each VOI configuration, the median value for all montages was , meaning that globally all montages exhibited good detectability performances. As expected, the DD montage showed overall best performance in terms of median and variance of AUC values. For focal VOIs [Fig. 8(a)], the and (respectively, and ) showed very similar performance in terms of median AUC values, suggesting our ability to reach a similar level of detectability with fewer sources and detectors. Results showed slightly more variability for the and montage than for the and montage. This interesting property suggests that after optimization of the sensor positions, the topographic maps were less dependent on the location of the VOI. For extended VOIs [Fig. 8(b)], and were showing very similar median AUC values, whereas the median value for the montage was slightly lower than the median value of the montage.
In simulations involving more than one VOI [Figs. 8(c) and 8(d)], AUC results obtained with the and the montages were exhibiting more variability than the regular montages. This was likely due to the fact that fewer sensors were available to accommodate each VOI in configuration involving multiple VOIs. For the montage, the large variance of AUC values observed in double VOI configurations [Fig. 8(c)] was likely due to the higher proportion of cases in which the branch-and-bound algorithm did not allocate positions for optodes over one VOI as discussed previously.
Figure 9 represents the distribution of the local energy scores defined from the GT area. For all VOI configurations, the montages showed slightly larger LE than the montages and similar variability. On the other hand, the montages showed clearly larger LE than the montages and the DD montages in all configurations. These results suggest that we can record localized high amplitude measurements over the activation area. The montage showed the largest LE scores likely due to the fact that it features a denser arrangement of SD pairs, thus allowing more flexibility.
Figure 10 represents the distribution of the ratio of visible voxels and the distribution of the mean overlap score within the visible voxels grouped by montage types. For all VOI configurations, the optimal montages showed larger ratios of visible voxels than the regular montages suggesting that a larger volume of the target could be detected by optimal montages. For single and double focal VOIs, of the volume was visible by the optimal and regular montages, whereas for extended VOIs, of the volume was visible. This decrease was mainly due to the fact that the additional voxels of the extended VOI were localized deeply in the brain and were, therefore, not visible by any SD configuration. The ratio of visible voxels with DD montages was lower for other montages mainly due to its lower density of SD measurements on the scalp. Similar to LE scores, the montage showed the larger mean overlap value (between 6 and 8), suggesting that the visible voxels are seen by several overlapping SD measurements.
Practicability in Epilepsy
Experimentally, even if enough sources and detectors were available from fNIRS system [e.g., the NIRx system (NIRx, USA) allows for 24 dual wavelength sources and 32 detectors], the installation of whole-head montages impose practical challenges. Indeed, the setup times required to attach a large number of optodes on the cap and to check the optical contact quality during the acquisition are limiting factors. In epilepsy, we aim at proposing routinely prolonged EEG/fNIRS recordings (several hours) to increase our chance to record spontaneous epileptic discharges or seizures. Therefore, we need practical and well-adapted solutions. A first solution would be to choose a subset of these regular whole-head arrangements. However, several challenges would still remain, such as choosing the location of the optode holders to select when targeting several VOIs, coping with the uncertainties of SD sensitivity that depend on the underlying anatomy, and taking into account the variability of the source/detector separation. DD montages cannot be used routinely because they cannot be integrated easily with a clinical arrangement of EEG electrodes.
The proposed method allows the experimenter to choose a limited number of optodes and provides automatically an adequate set of positions for optodes on the EEG/fNIRS cap. The algorithm may be applied easily to any EEG/fNIRS cap. For instance, the NIRx company (NIRx, USA) proposed a cap that can accommodate up to 256 electrodes or optodes following an extended 10/20 EEG coordinates system.69 Ishikawa et al. proposed a whole-head DD configuration70 with 128 optode holders and 64 EEG electrode positions interleaved between them.
Our proposed strategy requires some planning. Indeed, the segmentation of the subject anatomical MRI, the digitalization of the EEG/fNIRS cap, the MC simulations, and the montage optimization need to be performed before the acquisition. First, in the context of the presurgical evaluation of patients with epilepsy, collecting and segmenting anatomical MRIs is a typical procedure required by most multimodal imaging investigation (notably EEG/MEG source localization) and would not be specific to EEG/fNIRS. Second, the digitalization of the optode holder positions using the Brainsight neuronavigation system is a fast procedure ( for the isometric EEG/fNIRS cap). Finally, even if modern graphical processing units allow drastic reduction of the computation time, running the MC simulations was the most time-consuming procedure of our proposed method. Because it is not straightforward to predict without prior information on optical sensitivity that a given cortical area is likely to be measurable, we do believe, however, that modeling the light transport in complex head geometries is a necessary step to perform on each specific subject involved in a study in order to get more insights about the origins of the measured fNIRS signal.
In a nonclinical context, when studying a particular brain function among a group of subjects, one can consider using a template MRI and a template of possible fNIRS holders to design an optimal SD montage for the study. We assume that our approach would still provide relevant information when planning an EEG/fNIRS investigation. Careful evaluation of the added value of our approach when using a template MRI and more simplified model of light transportation would probably be sufficient in such a nonclinical context. We think that our algorithm could accurately provide optimal sets of optode positions in these conditions, but such an evaluation was, however, out of the scope of the present study.
Optimal Versus Regular Montages
In this study, we compared the montages computed by the optimization algorithm to regular whole-head montages with many sources and detectors (60 sources and 63 detectors for the montage, 92 sources and 156 detectors for the montage). Overall, optimal montages showed similar median AUC values and more variability of AUC distributions than our regular montages. However, all AUC values remained within an acceptable range (the median of the distribution was for multiple VOIs or 0.9 for single VOIs). For both optimal and whole-head regular montages, the distorted shape of some topographic maps resulted mainly from an interpolation artifact caused by the irregular spatial sampling of measurement points. This problem could be solved using different interpolation schemes. In the case of single VOI configurations (focal or extended), where only a subset of the SD pairs on whole-head montages are sensitive to the target VOI, our results suggest that the optimal montages have resolution characteristics similar to those of regular arrangements. In the case of multiple VOIs, we showed that even when using fewer sources and detectors for each target VOI, the estimated optimal montages were able to provide only minimal loss in detectability.
The optimal montages appeared more sensitive than regular montages. First, in some cases, the regular montages were not able to recover the full extent of the GT map (especially the montage) or were not able to generate topographic maps with equal contrast for multiple VOIs. This was mainly due to the fact that with regular montages the optode positions were not optimized; hence, some SD measurements were not sensitive to the target VOI even when located just above it.
Second, optimal montages showed larger LE values in the GT area than the regular montages. Even the DD montages that have good spatial resolution characteristics were exhibiting lower LE values. These results suggest that high-amplitude OD variations specific to some targeted brain areas are more likely to be captured with optimal montages. Because fNIRS signals are generally contaminated by different systemic physiological interferences coming from the extra cerebral layers,71 maximizing the sensitivity over cortical areas will tend to improve the SNR of the SD measurements. It means that less physiological interference should be present in the measured signals. Results in epilepsy mainly reported high SNR hemodynamic responses to seizure or generalized spike and wave,30,33,34 but very few studies reported responses to epileptic spikes (brief transient discharge not associated with clinical manifestations).31,32 It was demonstrated in these studies that statistical analysis was far from trivial, especially because of the interference from systemic fluctuations. Therefore, improving the SNR seems mandatory for spike analysis and actually we are more concerned about the sensitivity of our measurements than the spatial resolution.
Compared to fMRI, the spatial resolution of fNIRS is limited. fNIRT, however, can provide additional information because of its higher temporal resolution (20 to 100 Hz) and its unique ability to discriminate between and , which can also provide a measure of the relative change in cerebral blood volume under the assumption of constant hematocrit.72 Relative change in cerebral metabolic rate in oxygen and cerebral blood flow may be derived from these quantities in order to investigate possible impairment of neurovascular coupling at the time of epileptic discharges.34
In this study, we did not calculate the and topographic maps using the modified Beer–Lambert law mainly because experimentally fNIRT cannot recover accurately and , especially when dealing with focal activations.73 The MC approach, however, can provide subject-specific estimate of the differential path length and thus improve the quantifications of and .74,75 Note that the unknown partial volume errors62 may be estimated a priori from the VOI.
It was demonstrated in DOT that high-density montages with overlapping sets of measurements taken near the activated cortical region can greatly improve the spatial resolution and quantitative accuracy of the tomographic images.13,20,21,23,24 Interestingly, the optimal montages provide the greatest density of measurement in the targeted VOIs and the highest SNR. Therefore, we may conclude that the optimal montages are well adapted for optical tomography and allow us to solve this ill-posed inverse problem under better conditions. In order to verify this hypothesis, we performed a tomographic reconstruction from the noise-free simulated at 830 nm of one focal VOI using the , , , and montages. To do so, a Tikhonov regularization was used following the approach of White.24 The reconstruction was constrained within the GM. Figure 11 shows the reconstructed in one coronal slice for the four montages. The original VOI is represented in green in Fig. 11(a). Even if these results are very preliminary, we obtained, with significantly fewer sensors, a DOT accuracy similar to that obtained with the regular montages.
It was demonstrated that the maximum entropy on the mean (MEM) technique is sensitive to the spatial extent of the underlying EEG/MEG sources.45 In epilepsy, on the basis of the EEG/fMRI studies, we expect extended hemodynamic activity43,47 and we anticipate that the MEM framework will be appropriate for providing DOT reconstructions of epileptic hemodynamic activity. In future studies, we intend to evaluate the performance of MEM in the context of DOT when using optimal montages.
In the context of epilepsy research, a specific montage has to be designed specifically for each patient on an EEG/fNIRS cap in order to visualize image-functional or anatomical regions of interest. Therefore, given a set of available optode holder positions on an EEG/fNIRS cap, we proposed an original method to find the best SD montage automatically. The proposed validation framework included a simulation of focal and extended cortical activity and a comparison of validation metrics across different types of montages. We found that the topographic maps of optimal SD arrangements had spatial resolution properties comparable to those of regular SD arrangements but SNR characteristics that are better than those of SD regular arrangements. Moreover, the optimal montages provided more overlapping measurements, which is an important characteristic when performing DOT.
This work was supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant Program (C.G. and J.M.L.) and a pilot grant from the Quebec Bioimagery Network (QBIN). Electroencephalography/near-infrared spectroscopy equipment was acquired using grants from NSERC Research Tools and Instrumentation Program and the Canadian Foundation for Innovation (C.G.). C.G. and E.K. were also supported by a salary award from the Fonds de Recherche en Santé du Québec (FRQS). A.M. was supported by the Industrial Innovation Scholarships from the Fonds Québécois de la recherche sur la nature et les technologies (FRQNT), NSERC, and the Rogue Research company (Montréal, Canada) joint program. We would like to thank Professor Louis Collins and Dr. Vladimir Fonov for providing us with the the T1 magnetic resonance imaging segmentation and classification method.
Alexis Machado is a PhD student in the McGill University Department of Biomedical Engineering. He obtained a master’s degree in biomedical engineering from the École Polytechnique de Montréal, Canada. His research interests are in statistical signal processing and localization of epileptic spikes using diffuse optical imaging.
Odile Marcotte is an adjunct professor in the Department of Computer Science of the Université du Québec à Montréal and deputy director of the Centre de recherches mathématiques (Université de Montréal). She obtained a master’s degree in operations research from the École Polytechnique de Montréal and a PhD in operations research from Cornell University. Her research interests are in graph theory, combinatorial optimization, and integer programming.
Jean Marc Lina received an electrical engineering diploma from the National Polytechnic Institute of Grenoble, France, in 1982 and a PhD in theoretical physics from the Université de Montréal, Canada, in 1990. He is currently a professor in the electrical engineering department of the École de Technologie Supérieure and a member of the Centre de Recherches Mathématiques, Université de Montréal. His current research interests include wavelet analyses, inverse problems, entropic inference, and probabilistic graphical models applied to electromagnetic neuroimaging.
Eliane Kobayashi obtained her medical degree at UNICAMP, Brazil, where she also completed the neurology residency program, the epilepsy fellowship, and PhD studies. She came to the Montréal Neurological Institute (MNI) for postdoctoral studies. She is an assistant professor in the Department of Neurology and Neurosurgery at McGill University. She works as an epileptologist and neurophysiologist at the MNI Epilepsy Service and EEG Department. She received the Early Career Physician Scientist Award from the American Epilepsy Society.
Christophe Grova received an engineering and a MSc degrees in biomedical engineering (University of Technology of Compiègne, France), followed by a PhD in medical imaging (University of Rennes, France). Then, he came at the MNI as a postdoctoral fellow. Since 2008, he is assistant professor in Biomedical Engineering and Neurology and Neurosurgery Departments (McGill University). His research aims at investigating the integration of multimodal functional data to study brain mechanisms during epileptic discharges.