2.1.

Head Atlas and Meshing

To achieve a simulation setting which is close to a real measurement, we used the Montreal ICBM 2009a atlas, an unbiased nonlinear average of 152 anatomical MR images with $1{\mathrm{mm}}^3$ voxel size,^27,28 and corresponding tissue probability maps for cerebrospinal fluid, gray matter, and white matter.²⁹ In order to obtain a five-compartment model including scalp and skull, we additionally segmented the ICBM2009a images using mathematical morphological operations.³⁰ Based on this segmented brain atlas [see Fig. 1(a)], we used a masking and meshing software (Nirview)³¹ to create a 3-D tetrahedral mesh [Fig. 1(c)]. This mesh was used to calculate the photon transport, and thus provides the framework to simulate the cortical activation and to test the outcome of different reconstruction methods.

Fig. 1

(a) Segmented head atlas (ICBM 2009a, a nonlinear average of 152 MR images). From outer to inner layers: scalp, skull, cerebrospinal fluid, gray matter, white matter, (b) sketch of the optical fiber setup as used in the forward model (first nearest neighbor distance: 13 mm, (c) finite-element (FE) mesh of the left hemisphere with optical properties (${\mu }_a$), (d) example of the total sensitivity of determined from the unconstrained Jacobian $J$. A cross section of the sensitivity volume is superimposed on the corresponding layer of the head model, (e) total sensitivity of the spatially constrained Jacobian $\tilde{J}$: sensitivities for skull, scalp and CSF were set to zero.

2.2.

Forward Simulation and Spatial Constraints

Optical fiber positions on the boundary of the FE mesh were chosen according to the setting of a previous real-world cerebral DOT experiment conducted under resting conditions [Fig. 1(b)]. Due to the use of registration landmarks from the EEG 10–20 reference system³² and known source-detector distances, the coordinates for each fiber were known. To model light propagation, we used the Nirfast software toolbox,³³ a MATLAB®-based publicly available light modeling and reconstruction software. Nirfast applies the diffusion equation approximation, which is appropriate when scattering events dominate over absorption and the medium can be assumed to be an isotropic fluence field.

One challenge in DOT is the sensitivity of the measurement to signals coming from noncortical regions. The ${\mathrm{HbO}}_2$ specific wavelength is often contaminated with hemodynamic fluctuations from superficial veins in the scalp.³⁴ On the other hand, the decrease in absorption from the HbR sensitive wavelength is highly correlated to the BOLD response in functional magnetic resonance imaging.³⁵ For this simulation study, we use light model and data from the HbR sensitive wavelength of 760 nm. Optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ were assigned to each node of the FE mesh according to Strangman et al.³⁶

The result of the simulated light propagation is a sensitivity/Jacobian matrix $J$ with dimensions $M\times N$, where $M$ is the number of measurements (optical data channels) and $N$ is the number of nodes in the reconstruction volume. $J$ describes the logarithmic relationship between changes in measured boundary data ($\mathrm{\Delta }y$) that are caused by small changes of ${\mu }_a$ within the tissue for each channel-node combination, where

Since the reconstruction volume was not entirely covered with the optode set, and since DOT has only a limited penetration depth of 3 to 4 cm, we constrained $J$ in order to reduce the result space and thus reduce the “degree of ill-posedness.” One criterion for the exclusion of nodes was their affiliation to the noncortical tissue. Nodes belonging to scalp, skull, or cerebrospinal fluid were discarded. To exclude “weak” channels with very low sensitivity (e.g., due to large source-detector separation), we calculated the vector norm for all rows of $J$. Rows having a norm lower than 1% of the maximum value were discarded. The same procedure was performed for the mesh nodes, excluding nodes from the result space that had hardly been reached by any measuring channel. This step reduced the result space from 256 to 232 channels, and from 150,000 to 10,000 nodes. In the following, we refer to this reduced Jacobian as $\tilde{J}$ with the dimension of $\tilde{M}$ measuring channels and $\tilde{N}$ reconstruction nodes. Figure 1(d) depicts the total sensitivity of $J$ and Fig. 1(e) depicts $\tilde{J}$, which is calculated as the sum of the sensitivity over all measurement pairs for all used nodes within the head volume.

2.3.

Signal Generation and Noise Model

As an input signal, we modeled a hemodynamic response function (hrf) for absorption changes at 760 nm peaking 5 s after stimulus onset,³⁷ thereby mimicking a 400 s experiment with a stimulus duration of 20 s and an interstimulus interval of 20 s. This was necessary for testing the LCMV beamformer reconstruction method, which requires time-series data. Moreover, it allowed us to superimpose the artificial data with realistic noise of the same dimensionality obtained from the abovementioned resting-state recording.

Detector readings were generated as follows. A sparse matrix $A_{\mathrm{sim}}$ with the dimension of $\tilde{N}\times {\tilde{N}}_{\text{active}}$ was created, where ${\tilde{N}}_{\mathrm{c}}$ is the number of “activated” nodes. Each column of $A_{\mathrm{sim}}$ labels one node by setting $A_{\mathrm{sim}(l)}=1$ at a specific location $l$, while all other nodes are set to “0.” The locations for these “activated” nodes were randomly chosen, but due to restrictions of the reduced Jacobian $\tilde{J}$, all nodes were cortical. The specific sensitivity pattern $p$ in the activated node/nodes is defined by

We applied a realistic noise model for the purpose of testing different reconstruction algorithms under natural conditions. Most studies added white noise to the data to simulate the real measurements. In real life measurements, however, the noise is usually temporally and spatially correlated and is not normally distributed. We typically observe a higher increased noise level for larger source-detector separations than that for short distances. Second, the noise often has a high fraction of hemodynamic oscillations, which may interfere with the hemodynamic response and are sometimes hard to remove. Rather than applying a random noise term, we utilized data from a 10-min DOT experiment conducted under resting conditions as the noise model $\aleph$. For recording these resting-state data, we used a compact tomography imager that provides up to $32\text{sources}\times 32\text{detectors}$ (NIRScoutX Tomography Imager, NIRx Medizintechnik, Berlin, Germany). This allowed us to achieve realistic simulation data with characteristic features of the real measurements. The setup for that resting-state experiment was the same as that of the simulation setup, so that fiber distances and orientations were preserved.

We selected the rows of the noise matrix $\aleph$ according to the choice of channels for $\tilde{J}$, so that identical channels were used. Additionally, we took a subset of sampling/time points (columns) from $\aleph$, so that $y$ and $\aleph$ had the same dimensions. Finally, $\aleph$ and $y$ were normalized by their respective Frobenius norms in order to calibrate the artificial measurement and noise matrix. Given $y$, $\aleph$ and $s$, where $s$ is the signal level with a value between 0 and 1, the noisy simulated measurement $y_{\aleph }$ was constructed as

According to the real measurements, we low pass-filtered (first-order Butterworth) the generated data with a cut-off frequency of 0.4 Hz to remove the cardiac signals. In Fig. 2(a), we see detector readings from the resting measurement for large, medium, and short source-detector separations, and the dependence of the noise level on the fiber distances. Figure 2(b) depicts examples of the generated signal for two different measurement channels, each with a signal strength of 50% ($s=0.5$). Because of different locations and source-detector separations, the signal in the upper measurement is less dominated by noise compared to the second example in the lower graph.

Fig. 2

(a) Simulated DOT measurement with additive realistic noise recorded in resting condition using a compact tomography imager (NIRScoutX, NIRx, Medizintechnik, Germany). The noise level strongly depends on the source-detector separation, (b) two different measurement channels and generated signals with no noise (blue line) and with 50% noise (green line) added. The lower channel is noise dominated, since the signal generated is 100-fold smaller compared to the upper example. Due to different location and source-detector separations, noise has a different impact on the generated signal and may hamper the correct reconstruction.

2.4.

Image Reconstruction Methods

Since the number of measurements is much smaller than the number of reconstruction nodes, the linear system of Eq. (1) is heavily underdetermined, and a unique solution for $\mathrm{\Delta }{\mu }_a$ can only be obtained under constraints on the absorption coefficient distribution. In order to find a solution which is neuro-physiologically plausible, these constraints should always encode valid prior assumptions on the properties of $\mathrm{\Delta }{\mu }_a$. Various such assumptions have been proposed in the literature on an EEG/MEG inverse problem, which has a similar mathematical structure. In the following paragraphs, we introduce the approaches tested. As in previous parts of the paper, we omit the dependence of $\mathrm{\Delta }{\mu }_a$ and $\mathrm{\Delta }y$ on time. Thus, unless stated otherwise, a separate reconstruction is performed for each time point (i.e., difference measurement).

2.5.

Reconstruction Quality Criterion: Earth Mover’s Distance

Each of the image reconstruction procedures resulted in a matrix with the dimension $\tilde{N}\times {\text{samples}}^{\mathrm{hrf}}$. To estimate the quality of the result, we calculated a general linear model in a sense of a linear regression for all reconstructed time courses $x_{1,\dots ,\tilde{N}}$ with hrf as the regressor. Thus, for each voxel of the reconstruction volume, a $t$-value was derived. All negative $t$-values and those with a $t$-value smaller than 20% of the maximum $t$-value were eliminated.

As a measure of overall reconstruction quality, we applied the Earth Mover’s Distance (EMD⁴²) to the reconstruction results ($t$-values) of all methods. The EMD calculates the minimal amount of energy that must be spent to transform one distribution into the other. Given the known locations ($xyz$-coordinates) of the mesh nodes in a 3-D space, the EMD uses the Euclidean distance between all nodes as a ground metric to calculate the minimum costs of transforming the normalized histogram of $t$-values into the normalized histogram of the simulated activations. Figure 3 shows an impression of a good reconstruction result with a low EMD [Figs. 3(b) and 3(c)] and a poor result [Fig. 3(d)] based on one simulated activation [Fig. 3(a)]. The advantage of the EMD is its ability to compare the overall distribution of the 3-D volumes. Unfortunately, solely looking at the EMD value gives no hint as to whether the result is blurry and/or dislocated. To gain additional information about the reconstruction quality in terms of the malpositioning of the activation, we additionally calculated the Euclidean distance between the simulated target and the maximum value of the result for the cases where only one spot was activated.

Fig. 3

Example for image reconstruction using tSVD and with different numbers of singular values used for inversion of $\tilde{J}$. (a) Simulated target activation, (b) result using 30 used singular values for reconstruction ($\mathrm{EMD}=12.6$, best possible result), (c) result using (cross-validated) 60 singular values ($\mathrm{EMD}=15.1$), (d) result from reconstruction with 160 singular values ($\mathrm{EMD}=57.3$).

2.6.

Automatic Determination of the Regularization Parameter Using Cross-Validation

Distributed inverses, such as ${\ell }_2$MNE, ${\ell }_1$MNE, ${\ell }_0$MNE, tSVD, wMNE, and S-FLEX, directly estimate the source activity $\mathrm{\Delta }{\tilde{\mu }}_a$. This means that for an $\tilde{M}\times T$ sensor time series, $\tilde{N}\times T$ parameters have to be estimated, where $\tilde{N}\gg \tilde{M}$. Under these circumstances, regularization is necessary (as outlined above), and the choice of the regularization parameter crucially affects whether the fitted model is too complex (overfitting the data), too simple (not explaining the relevant aspects of the data), or “just right.”

Beamformers, on the other hand, are characterized by a low number of parameters. Therefore, the estimation is typically very stable. The LCMV beamformer in Eq. (19), for example, solves $\tilde{N}$ optimization problems (one for each voxel), each of which is concerned with the estimation of only a single $\tilde{M}$-dimensional filter ${\tilde{v}}_q$ based on the covariance matrix of a $\tilde{M}\times T$ dataset $\mathrm{\Delta }{y}_{\aleph }$, where $T$ is the number of samples, and typically $T\gg \tilde{M}$.

The parameter $\lambda$ of regularized models drives the estimated brain activation ($\mathrm{\Delta }{\tilde{\mu }}_a$) away from the solution that explains the best measurement to a solution with “simpler” structure. As such, $\lambda$ critically affects the shape of the chosen solution and the reconstruction accuracy. Therefore, choosing the “right” amount of regularization is very important. This choice should not be based on visual inspection or other subjective measures in order not to bias the later neurophysiological interpretation of the results. Rather, an automatic selection criterion is required.

One way of assessing the quality of a regularized model is to measure how well it explains the unseen data which have not been used for estimating the model parameters. This can be done using CV. In $k$-fold CV, the data are split into $k$ chunks. The model is fitted on $k-1$ chunks and evaluated on the remaining “test” chunk. This procedure is repeated for each choice of the regularization parameter and for each choice of the test chunk. The parameter that best explains the test data on average is selected is and used for training a final model based on the entire data available.

In the distributed inverse source reconstruction, data folds are created by dividing the measurement channels into $k$ sets, and the performance criterion to be estimated is the squared loss at the “test” channels, i.e., ${\Vert {\tilde{J}}^{\text{test}}{\tilde{\mu }}_a-\mathrm{\Delta }{y}_{\aleph }^{\text{test}}\Vert }_2^2$, where ${\tilde{J}}^{\text{test}}$ and $\mathrm{\Delta }{y}_{\aleph }^{\text{test}}$ are the parts of $\tilde{J}$ and $\mathrm{\Delta }{y}_{\aleph }$ belonging to the test channels.

For inverse methods estimating the brain activations as linear combinations of the data using some pseudoinverse ${\tilde{J}}_{\lambda }^{\#}$ (such as MNE, wMNE, and tSVD), an approximation to leave-one-out CV (i.e., $k$-fold CV with $k=M$) can be carried out in the closed form. The so-called generalized CV score $g(\lambda )$ is given by

One goal of this work is to show how the reconstruction quality alters when different regularization values are used for reconstruction. Methods that directly estimate $\mathrm{\Delta }{\mu }_a$ are highly dependent on the choice of this parameter. To visualize this relationship, we first generated one target, then added 50% noise to the artificial measurement matrix, and finally reconstructed this specific activation using a wide range of values for $\lambda$. For every instance of this reconstruction result, the EMD was calculated. This procedure was repeated 50 times for ${\ell }_2$MNE and wMNE. To test the same for tSVD, we proceed in the same manner except that we increased the number of singular values used for reconstruction, starting with the 10 highest and ending with using all ($m=231$).

3.1.

Reconstruction Quality Highly Depends on the Choice of the Regularization Parameter: an Almost Optimum Choice Can be Made without a User Bias

To visualize the impact of the chosen value for regularization, Fig. 3 depicts an example of reconstruction for tSVD, where the activation was recovered using different numbers of singular values for the inversion of $\tilde{J}$. Figure 3(b) shows the best possible reconstruction result with the lowest EMD for this simulation [Fig. 3(a)]. The result that was achieved with the cross-validated number of singular values is shown in Fig. 3(c). Both parameters resolve the activation reasonably well with a correct location and little blur. The result obtained with 160 used singular values [Fig. 3(d)] leads to overfitting, which is evident from the high number of phantom activations.

Figure 4 depicts multiple graphs, each representing one of the distributed reconstruction methods used. The red solid line shows the mean EMD for 50 different simulations and a wide range of values for $\lambda$ (increasing number of singular values $m$ for tSVD, respectively). The red transparent area represents the standard error of the mean and the blue area represents the standard deviation. In all quality plots, we clearly see how EMD changes with different regularization parameters. We find a high EMD when very small or very high regularization values are chosen, rendering data that are either over or under fitted. Between them, we find a global minimum, which is indicated by the red dot representing the best possible EMD. Assuming that the location of this minimum would be known prior to reconstruction, this $\lambda$ ($m$, respectively) would be the first choice for parameter selection. However, in real-world experiments, the true location and extent of this activation is unknown and such a plot is not available. To overcome this challenge, this optimum is approximated by the CV as described in the section above. The blue dots in each subplot indicate the mean value for $\lambda$ ($m$, respectively), estimated using the CV and the respective mean EMD. In all three methods, the cross-validated value leads to results that are comparable in quality to the best possible result. The slight mismatch between the best possible and cross-validated results may be caused by the limited amount of data available.

Fig. 4

Depiction of the relation between regularization and reconstruction quality of three distributed reconstruction methods (noisy data, one activated spot). (a) Result for ${\ell }_2$ MNE. In each simulation run, the activation was reconstructed using 100 different regularization parameters. The red line depicts the average EMD for 50 simulation runs. The geometric mean of the best possible regularization value (red dot) and the same for the automatically detected (cross-validated) (blue dot) and their respective mean EMD. (b) Reconstruction quality for tSVD using an increasing number of singular values for reconstruction, (c) result for wMNE.

Please note that since ${\ell }_1$MNE, ${\ell }_0$MNE, and S-FLEX cannot be solved in the closed form and rely on numerical optimization, the calculation time for such a large number of variations was unreasonably high. Therefore, the results for these methods are not shown here. In practice, we choose the regularization strengths for these methods indirectly by selecting $\lambda$ such that the data are explained to the same extent as it was explained by wMNE using a cross-validated $\lambda$. The LCMV beamformer is also omitted here, since it does not depend on the choice of a regularization parameter in the same way as do the other methods, as mentioned above.

With respect to the reconstruction quality concerning different amplitudes of the simulated target, we calculated additional simulation, testing two more aspects. First, we reconstructed a target on a fixed location with two different amplitudes and a fixed regularization parameter (optimized for one of the simulations). Then, we reconstructed a target in a fixed location with different amplitudes and a variable regularization parameter. In both cases, the difference in amplitude heights in the most “active” voxels reflected the simulated difference. The reconstruction quality was almost identical.

3.2.

Linearly Constrained Minimum Variance Beamforming Resolves Single Activation Spots Best

The second focus of this work is on benchmarking source reconstruction methods, among which are frequently used methods such as the recently proposed sparse algorithms and EEG-source localization methods. These methods are introduced below in the context of cerebral DOT with its semi-infinite geometry. Figure 5 gives an impression of the simulation and the reconstruction result for a single spot activation in a single case with the seven reviewed methods. For visualization, we show the transverse cross sections covering the area of the simulated activation.

Fig. 5

Exemplary reconstruction images for a single spot activation. (a) Transversal slices of the reconstruction volume with the simulated activation in column 6. The other columns depict transverse cross sections adjacent to the central layer ($z$ direction, slice depth: 1 mm). (b) Reconstruction result for a 0% noise level. Each row represents the result from one particular reconstruction method. The number in the right column indicates the Earth Mover’s Distance (EMD) for this specific example. ($c$) 50% noise added to the data.

The arrow in Fig. 5(a) indicates the node that was set “active.” Rows 1 to 7 in Fig. 5(b) show the reconstructed images for all tested methods in a noise-free simulation. Within each row, the EMD between the simulation and the result is pointed out in the last column. Figure 5(c) shows the same simulation but with 50% noise in the data.

For ${\ell }_2$MNE, tSVD, and wMNE, we find a relatively good localization of the peak activation with slight blurring in the noise-free simulation. This blurring increases when noise is added to the data. Compared to ${\ell }_2$MNE, wMNE shows an increased sparsity and a lower EMD. S-FLEX and ${\ell }_1$MNE show small positioning errors in the noisy case and a focal result. In both noise levels, we find the ideal results for LCMV, with no displacement and a high focality. All the three latter methods appear to be rather insensitive to the applied noise level. ${\ell }_0$MNE performs well in the noise-free case, but fails when noise is added to the data.

For an overall comparison of all methods, the average EMD of 100 simulations with one activated spot and four different noise levels (0%, 25%, 50%, and 75%) can be found in Fig. 6(a) The respective mean Euclidean distance between simulation and maximum value of the reconstruction result can be found in Fig. 6(b).

Fig. 6

(a) Overall EMD statistics for single spot activation, four applied noise levels and all seven reconstruction methods, (b) averaged Euclidean distance between simulated target and maximum value of the result in mm for all methods and noise levels. Black bar indicates the standard error of mean.

Similar to the single case, we find the best reconstruction at every noise level when LCMV is used. In almost all simulations, the beamformer achieves a correct positioning with minimal blurring even at the highest noise level. S-FLEX and ${\ell }_1$MNE perform well and recover sparse results; however, their results are dislocated by a few millimeters. Interestingly, S-FLEX and ${\ell }_1$MNE do not achieve their best EMD scores at the lowest noise levels with high signal levels [see Eq. (4)]. This may be due to the fact that, for efficiency reasons, the optimization for both methods is stopped after the data have been fitted with a goodness-of-fit of $\mathrm{gof}=0.95$, where $\mathrm{gof}=1-{\Vert \tilde{J}\mathrm{\Delta }{\tilde{\mu }}_a-\mathrm{\Delta }{y}_{\aleph }\Vert }_2/{\Vert \mathrm{\Delta }{y}_{\aleph }\Vert }_2$ The data may be insufficiently fitted for very low noise levels.

TSVD, ${\ell }_2$MNE, and wMNE show a clear dependence on the noise level: with higher noise, the EMD increases. This can be especially observed in tSVD. However, although reaching a high EMD, tSVD still shows only a small positioning error (Euclidean distance) between the peak value of the reconstruction and the simulation [Fig. 6(b)]; within the highest noise level, the average Euclidean distance between the result and simulation is 8.3 mm (${\ell }_2$MNE: 15, wMNE: 11, LCMV: 0.2, S-FLEX: 10.1, and ${\ell }_1$MNE: 8.8 mm).

This implies that the main reason for a high EMD is a higher blur level rather than malpositioning; this blur could possibly be reduced by thresholding the result. The highest sensitivity to noise is found in ${\ell }_0$MNE: beginning with low noise levels, the EMD and the positioning error dramatically increase.

3.3.

Minimum ${\ell }_1$-Norm Achieves Best Result When Two Spots Are Active

When investigating a relatively small area of the brain, there is often only one spot of activation within the probe. However, there are approaches where larger areas or even the whole head are scanned. When the medium is larger, the possibility of including two or more areas with simultaneously fluctuating rhythms caused by a synchronic hemodynamic answer rises. We, therefore, simulated two additional areas with perfectly correlated activity in the brain.

Recovering two (or more) activation foci in an algorithm is a challenge. TSVD, ${\ell }_2$MNE, and wMNE show no significant differences in their EMD, which is attributable to the generally increased level of blur. That makes it harder to distinguish the quality using the EMD method. However, when looking at single cases with visualized reconstruction results (Fig. 7), we can see that all methods except the beamformer are capable of recovering both activations. Since ${\ell }_1$MNE can reconstruct the sparser activation patches more clearly than the other methods, its performance is better, although again some slight positioning errors do occur. At lower noise levels, S-FLEX yields results comparable to those of ${\ell }_1$MNE, but their quality decreases at the highest noise levels. ${\ell }_0$MNE can almost perfectly recover both targets in a noise-free dataset, but fails again when noise is added. Due to reduced blur, wMNE shows a slight but not significant advantage over ${\ell }_2$MNE, and with increased noise levels it also has a slight advantage over tSVD. Finally, it is obvious that the LCMV beamformer cannot resolve correlated activity at different brain sites, and, therefore, shows a greatly decreased performance. For a comparison see Fig. 8.

Fig. 7

Exemplary reconstruction images for two activations. (a) Simulated activation: two nodes in different locations were defined as “active” (indicated by the arrow), (b) reconstruction results for noise-free data and from seven different reconstruction algorithms, (c) results for noisy data (50%). Columns represent transverse cross sections of the reconstruction volume ($z$-direction, slice depth: 1 mm).

Fig. 8

Overall EMD statistics ($n=100$) for all seven methods and four different noise levels and two activated spots. Black bar indicates the standard error of mean.

3.4.

Test on Experimental Data

It is always of interest to see how algorithms work with real data. However, it is difficult to estimate or compare the reconstruction quality with no objective reference. To get an impression of how the different algorithms work with real data, we added reconstruction results for a choice of reconstruction methods (LCMV, ${\ell }_2$MNE, tSVD, and wMNE) for a finger tapping task (right hand tapping for 20 s followed by 20 s rest, 10-min duration). Activated areas were identified using a general linear model. In Fig. 9, we show the lateral view on the left hemisphere with colored areas indicating voxels with a significant ($t$-values $\le -4$) hemodynamic answer in the HbR time courses.

Fig. 9

Reconstruction result for a choice of methods [(a) LCMV, (b) ${\ell }_2$MNE, (c) tSVD, (d) wMNE] and a finger tapping task of the right hand side. Lateral view on the left hemisphere. Colored areas indicate the voxels with a significant hemodynamic response in the HbR time courses due to the finger tapping (estimated with a GLM, $t$-values $\le -4$).