Short-separation regression incorporated diffuse optical tomography image reconstruction modeling for high-density functional near-infrared spectroscopy

Yuanyuan Gao; De’Ja Rogers; Alexander von Lühmann; Antonio Ortega-Martinez; David A. Boas; Meryem A. Yücel

doi:10.1117/1.NPh.10.2.025007

23 May 2023 Short-separation regression incorporated diffuse optical tomography image reconstruction modeling for high-density functional near-infrared spectroscopy

Yuanyuan Gao, De’Ja Rogers, Alexander von Lühmann, Antonio Ortega-Martinez, David A. Boas, Meryem A. Yücel

Author Affiliations +

Neurophotonics, Vol. 10, Issue 2, 025007 (May 2023). https://doi.org/10.1117/1.NPh.10.2.025007

Abstract

Significance

Short-separation (SS) regression and diffuse optical tomography (DOT) image reconstruction, two widely adopted methods in functional near-infrared spectroscopy (fNIRS), were demonstrated to individually facilitate the separation of brain activation and physiological signals, with further improvement using both sequentially. We hypothesized that doing both simultaneously would further improve the performance.

Aim

Motivated by the success of these two approaches, we propose a method, SS-DOT, which applies SS and DOT simultaneously.

Approach

The method, which employs spatial and temporal basis functions to represent the hemoglobin concentration changes, enables us to incorporate SS regressors into the time series DOT model. To benchmark the performance of the SS-DOT model against conventional sequential models, we use fNIRS resting state data augmented with synthetic brain response as well as data acquired during a ball squeezing task. The conventional sequential models comprise performing SS regression and DOT.

Results

The results show that the SS-DOT model improves the image quality by increasing the contrast-to-background ratio by a threefold improvement. The benefits are marginal at small brain activation.

Conclusions

The SS-DOT model improves the fNIRS image reconstruction quality.

1. Introduction

Continuous-wave functional near-infrared spectroscopy (fNIRS) is a noninvasive neuroimaging technique that measures cerebral hemoglobin concentration changes using an array of light sources and detectors placed on the surface of the head.¹^–⁵ The detected light reflects the hemodynamic changes not only in the cortex but also in the superficial tissues through which the light propagates, i.e., scalp, skull, and cerebrospinal fluid (CSF).⁶^–⁹ Moreover, measurements are more sensitive to extracerebral tissues that are closer to the detectors and much less sensitive to the brain cortex, where the brain functional activity takes place.¹⁰ Thus accurately modeling the cerebral hemodynamic response function (HRF) and distinguishing it from the extracerebral signals is essential to robustly estimating brain activation.

Various methods to filter the systemic interference from the fNIRS measurement have been described in the literature. Low-pass filtering is widely used as it is highly effective at removing high-frequency components in the signal such as cardiac oscillations.⁴^,¹¹ However, it tends to overcorrect the signal when the brain signals have a significant overlap in the frequency spectrum with other physiological variations, such as respiration or Mayer waves.¹² As such, more powerful methods for global noise reduction have been developed. For example, auxiliary physiological measurements are regressed out from fNIRS measurements, such as heart beat,¹³^–¹⁵ blood pressure and heart rate,¹⁵^–¹⁷ scalp hemodynamic changes obtained by fMRI,¹⁸ and another source–detector measurement performed over a nonactivated region of the brain.⁴ Data-driven methods were also used to reduce global noise, such as principal component analysis¹⁹^,²⁰ and independent component analysis.²¹

Beyond the previously mentioned methods, two fundamental approaches have been suggested: (1) inclusion of short-separation (SS) measurements as a systemic physiological regressor and (2) high-density diffuse optical tomography (DOT) to spatially resolve brain activity from physiological contamination in the scalp. The SS measurement approach²²^–²⁷ utilizes the fact that long-separation channel ( $\sim 30 mm$ ) signals sample both cerebral and extracerebral tissues and SS channel ( $\sim 8 mm$ ) signals mostly sample the extracerebral tissues. The signal measured with the SS channels, representing mostly the extracerebral hemodynamics, can be regressed out from the signal measured at the long-separation channels, improving the estimation of cerebral hemodynamics. Moreover, some of the systemic physiological interference arising from cardiac activity, respiration, and other homeostatic processes that appear both in cerebral and extracerebral layers can also be captured by the SS measurements and can be regressed out from the long separation measurements. Two widely used methods for the estimation of the HRF are the linear least square estimator/general linear model (GLM) approach²²^,²⁷^–²⁹ and the adaptive filtering technique that updates the coefficients at each time point.²³^,²⁷^,³⁰

DOT reconstructs three-dimensional (3D) maps of the hemoglobin concentration changes.³^,⁷^,⁹^,³¹ It is usually used on datasets collected with a high-density probe consisting of a dense array of sources and detectors that provide overlapping measurements typically with channels consisting of multiple distances.⁹ The overlapping measurement can increase both the spatial resolution and the depth resolution. The disentanglement of extracerebral and cerebral hemodynamics is realized by the varying depth sensitivity of the measurements with different source–detector distances, as long-distance channels have more relative sensitivity to the deeper cerebral layer and the short-distance channels are more sensitive to the extracerebral layers. In contrast to the SS regression approach, which is typically applied in channel space, DOT image reconstruction results in an image of the hemoglobin concentration changes. DOT image reconstruction technique has evolved in the past 20 years,⁹^,³²^–³⁴ decreasing the localization error and improving the effective resolution of the images.³⁵ To construct the DOT forward model, finite-element³⁶^–³⁸ or Monte Carlo methods³⁹^–⁴² can be employed. To solve the inverse problem, the regularization priors have also been optimized for DOT, such as depth sensitivity regularization,⁴³ singular value decomposition regularization,⁴⁴ and regularization methods adapted from EEG source reconstruction like linear constrained minimum variance beamforming.⁴⁵

One can either estimate the HRF in channel space and reconstruct a brain map using the resultant HRF amplitudes⁴⁶^,⁴⁷ or first reconstruct an image time series from the fNIRS channel data and then estimate the HRF in image space.⁹^,⁴⁸ Gregg et al.²⁶ tested the performance of implementing SS and DOT approaches by applying them sequentially. The contrast-to-noise ratio was increased, and the image variance across blocks was decreased with either one of the methods, but the best performance was achieved when sequentially combining the two. Simultaneously reconstructing the image and estimating the HRF have also been suggested in the literature.¹⁷^,⁴⁹^,⁵⁰ This approach uses temporal and spatial basis functions to simultaneously model HRF over time and space. The model is flexible and can incorporate auxiliary inputs as regressors, such as arterial blood pressure, respiration, or motion. Motivated by these works, we proposed merging the SS and DOT models into one model and solving it simultaneously. We hypothesize that it would further improve the contrast-to-background ratio (CBR), $t$ -value, and lateralization index.

In this work, we establish a DOT model that incorporates SS regression (SS-DOT) to reconstruct brain activation maps. We design a high-density probe, by which the spatial and temporal basis functions enable simultaneous spatial and temporal reconstructions of the brain activation. We employ Monte Carlo simulation to derive the sensitivity matrix. We then construct the design matrix to comprise brain activation, extracerebral stimulus activation, SS regressors, and drift regressors to set up the forward and inverse models. Using data from simulation and experimental studies, we test SS-DOT against the conventional sequential models, i.e., performing SS regression first and then the DOT image reconstruction. We quantify the performance of each model using the metrics of CBR, $t$ -value, and lateralization index.

2. Methods

2.1.

High-Density Probe Layout

The high-density fNIRS probe used in this study was designed using a hexagonal pattern [Fig. 1(a)], with 7 sources and 16 detectors forming 50 channels. The SS detectors are placed 19 mm apart from the sources, and the long-separation detectors are placed 32.9 mm apart from the sources, covering a $76 \times 66 {mm}^{2}$ area. This probe was placed on the motor region of the head, with only one patch on the left hemisphere in the simulation study and two patches on both hemispheres in the experimental study [Fig. 1(b)]. Following our “NinjaCap” approach, we 3D-printed the caps designed in AtlasViewer⁴⁷ to ensure the proper positioning of the optodes on the head.⁵¹ We also measured the Cz position of the head to make sure the cap’s Cz marker and the head’s true Cz position were accurately aligned.

Fig. 1

Probe design and illustration of the spatial kernels. (a) The layout of one patch of the high-density probe; red numbers are sources, and blue numbers are detectors. Green lines are channels. (b) The position of the probe patch on the head. (c) The sensitivity profile in the brain (left) and scalp (right) region on a log 10 scale. (d) The illustration of the 3D spatial Gaussian kernels. (e) The center positions of the spatial kernels of the brain (left) and the scalp (right) regions.

2.2.

Spatial and Temporal Bases

As described in Ref. 49, to represent the underlying physiological changes in both the brain and scalp regions, we assumed that changes in oxy-(HbO) and deoxyhemoglobin (HbR) concentration can be described as linear combinations of a set of spatial and temporal basis functions designated to capture the functional and physiological hemodynamic fluctuations. One of the benefits of implementing the spatial basis is to reduce the degrees of freedom of the model to make the model more computationally tractable.

We employed a set of overlapping 3D Gaussian kernels ( $g$ ) to represent the individual functional hemodynamic changes at each vertex of the brain and scalp surface [Fig. 1(d)]:

Eq. (1)

g = \exp (\frac{- {‖ x - x_{0} ‖}^{2}}{2 σ_{g}^{2}}),

where

x_{0}

represents the center point location of the Gaussian kernel,

x

represents the location of the element vertex, and

σ_{g}

is the Gaussian standard deviation, which controls the width of the Gaussian kernel. The spatial bases were placed only in the regions of the head with a high-sensitivity profile (i.e., sensitivity value

> 0.01

) to reduce computational expense. The kernels were placed 5 and 20 mm apart from each other for the cortex and the scalp, respectively [Fig. 1(e)]. The Gaussian kernel standard deviation

σ_{g}

was set equal to the kernel distance. The brain and scalp spatial basis functions were restricted to the element vertices on the brain and scalp, respectively. The spatial basis sets for HbO and HbR changes in brain and scalp are denoted as

G_{brain, HbO}

,

G_{brain, HbR}

,

G_{non - brain, HbO}

, and

G_{non - brain, HbR}

, respectively.

To model the temporal dynamics of the evoked functional changes within each of these Gaussian bases, we used a linear combination of a two-parameter modified Gamma function, as given by the following equation:

Eq. (2)

t (t) = \frac{{(t - τ)}^{2}}{σ^{2} e^{- 1}} \exp (\frac{- {(t - τ)}^{2}}{σ^{2}}),

where

τ

is the delay time of the response and

σ

is the parameter that controls the temporal width of the response function. This basis support is similar to the temporal basis often used in GLM analyses.⁵²^,⁵³ The overall temporal model of the functional component of the hemodynamic signals is expressed as the convolution of the experimental stimulus timing (

U

) and the functional impulse response function:

Eq. (3)

T = t * U,

where

U

is a binary vector describing the experimental paradigm (i.e., the timing of stimulus presentation) and spans the temporal duration of the experiment and * is the convolution operator. The dimensions of the

T

matrix are the number of temporal basis functions by the number of measurement time points. Here we used the same temporal basis function for the scalp region [Eq. (2)] to represent the physiological responses to the stimulus in the scalp region. The temporal basis sets for HbO and HbR changes in brain and scalp are denoted as

T_{brain, HbO}

,

T_{brain, HbR}

,

T_{non - brain, HbO}

, and

T_{non - brain, HbR}

, respectively.

Finally, all spatial–temporal bases are combined into a single Kronecker operator ( $GT$ ) with a total of four unique basis groups. This matrix is formed by the Kronecker product of the individual spatial ( $G$ ) and temporal ( $T$ ) basis functions for each of the eight groups mentioned above. The matrix $GT$ is given by the following equation:

Eq. (4)

GT = [\begin{matrix} G_{brain, HbO} \otimes T_{brain, HbO} & G_{non - brain, HbO} \otimes T_{non - brain, HbO} \\ G_{brain, HbR} \otimes T_{brain, HbR} & G_{non - brain, HbR} \otimes T_{non - brain, HbR} \end{matrix}],

which represents the hemoglobin changes in each vertex and each time point for each condition of stimulus. Here

\otimes

represents the Kronecker product.

2.3.

Sensitivity Matrix

We employed Monte Carlo simulations to calculate the sensitivity matrix $A_{λ}$ using the MCXLAB toolbox⁵⁴ for each wavelength, in this case, for two wavelengths $λ_{1}$ and $λ_{2}$ . The $A_{λ}$ matrix has the dimension of the number of channels by the number of voxels or vertices, and it maps the changes of the absorption coefficient values in voxels to the change in optical density values in the measurement channels. The Colin head model was used in this Monte Carlo simulation, and $10^{7}$ photons were simulated from each optode. The scattering anisotropy $g$ does not change the result in the diffusion equation provided that the value of the reduced scattering coefficient is kept constant. The Monte Carlo simulation is more efficient for smaller values of $g$ , so we used $g = 0.01$ .³³ We used a refractive index of 1 for all tissues.³³ Colin27 was used as our head model.⁴⁷^,⁵⁵ The Colin head model⁵⁵ was segmented into five types of tissues, scalp, skull, CSF, gray matter, and white matter using SPM software.⁵⁶ The absorption coefficients $μ_{a}$ and reduced scattering coefficients $μ_{s}^{'}$ used for simulations are listed in Table 1.⁵⁷^–⁶⁰

Table 1

Absorption coefficients μa and reduced scattering coefficients μs′ for different tissue layers.

Tissues	μa (mm−1)		μs′ (mm−1)
Tissues	760 nm	850 nm	760 nm	850 nm
Scalp	0.0177	0.0190	0.73	0.64
Skull	0.0125	0.0139	0.93	0.84
CSF	0.0021	0.0040	0.30	0.30
Gray matter	0.0195	0.0192	1.18	0.67
White matter	0.0195	0.0208	1.18	1.01

The foby

Eq. (5)

A = (\begin{matrix} ϵ_{λ_{1}}^{H b O} A_{λ_{1}} & ϵ_{λ_{1}}^{H b R} A_{λ_{1}} \\ ϵ_{λ_{2}}^{H b O} A_{λ_{2}} & ϵ_{λ_{2}}^{H b R} A_{λ_{2}} \end{matrix}),

where

ϵ

is the extinction coefficient for the chromophores (HbO and HbR) at each wavelength

λ_{i}

.⁶¹

2.4.

Forward Model

We first converted the light intensity into the changes in optical density values by

Eq. (6)

Δ OD (t) = - \log (ϕ (t) / ϕ_{0}),

where

Δ OD

is the change in the optical density, the logarithm is a common log, and

ϕ_{0}

is the average detected photon fluence. Before the conversion, we prune the noisy channels with a signal-to-noise ratio (SNR) lower than 5 using the hmrR_PruneChannels function in Homer3. The SNR is defined as the mean divided by the standard deviation of the signal. After the conversion, motion artifacts are corrected by the SplineSG method using the hmrR_MotionCorrectSplineSG function in Homer3, with the interpolation parameter set to 0.99.⁶²

The measurement from the channels was modeled as a linear combination of four components: (1) stimulus-derived cerebral activity, (2) stimulus-derived extracerebral activity, (3) global systemic physiology, and (4) drift in the measurement itself. The first two parts are described as a linear sum of the spatial–temporal basis function $GT$ weighted by a vector of unknown coefficients denoted by $b_{f}$ :

Eq. (7)

[\begin{matrix} Δ HbO \\ Δ HbR \end{matrix}] = GT \cdot b_{f} .

The optical density changes are obtained by

Eq. (8)

Δ {OD}_{f} = A \cdot GT \cdot b_{f},

where the subscript

f

represents the functional changes.

Global systemic physiology is captured by the measurement in the SS channels denoted by $Δ {OD}_{SS}$ , and the signal drift is represented by drift regressors denoted by ${OD}_{drift}$ . We denote the measurement from channels in the high-density probe (including the long 32.9 mm channels and the shorter 19 mm channels) as $Δ OD$ . Consequently, the forward model is written as

Eq. (9)

Δ OD = [Δ {OD}_{SS}, {OD}_{drift}, A \cdot G T] \cdot b,

where

b

are the unknown coefficients for SS regressors, drift regressors, and the functional changes.

This final model is written as a single matrix equation:

Δ OD = [Δ {OD}_{SS}, {OD}_{drift}, H] \cdot b,

Eq. (10)

Y = O \cdot b,

where

Y

is the concatenated set of measurement channels for all time points, i.e.,

Δ OD

;

O

is

[Δ {OD}_{SS}, {OD}_{drift}, H]

; and

H

is

A \cdot G T

.

2.5.

Inverse Model and Regularization

First, we employ spatial regularization,⁹ which is realized by replacing $A$ in Eq. (9) by $\hat{A}$ . $\hat{A}$ is the $A$ matrix spatially regularized by the diagonal matrix $L$ :

Eq. (11)

\hat{A} = {A L}^{- 1},

Eq. (12)

diag (L) = \sqrt[2]{diag (A^{T} A) + η} .

To estimate the coefficients of the spatial–temporal basis sets, the linear model Eq. (10) is solved using the Tikhonov regularization scheme (Moore–Penrose generalized inverse), with spatial regularization:⁹

Eq. (13)

b = L^{- 1} (O^{T} Y) / (O^{T} O + α I_{f}),

where

α

is the regularization parameter and

I_{f}

is a diagonal matrix that has diagonal values of ones for the regressors of

A \cdot G T

and zeros for the regressors of

Δ {OD}_{SS}

and

{OD}_{drift}

as we only regularized the functional changes to balance the high spatial frequency noise that would otherwise appear in the image because of the ill-conditioned

A \cdot G T

matrix. The Tikhonov parameter

α

adjusts the balance between high spatial frequency noise and image smoothness and was set to 0.01 times the maximum diagonal value of

O^{T} O

. The spatially variant regularization parameter

η

was set to 0.01 times the maximum diagonal value of

{\hat{A}}^{T} \hat{A}

.⁶⁰^,⁶³ The GPU acceleration method for this model is presented in the Supplementary Material.

2.6.

Sequential Models

The sequential model refers to a model that first performs SS regression in the channel-wise space and then carries out the image reconstruction.⁸ To do the regression of SS in the channel-wise space, we use the modified Beer–Lambert law to convert the measured OD time series data from each channel to Hb signals:

Eq. (14)

(\begin{matrix} {Δ OD}_{λ_{1}} \\ {Δ OD}_{λ_{2}} \end{matrix}) = (\begin{matrix} ϵ_{λ_{1}}^{HbR} ppf & ϵ_{λ_{1}}^{HbO} ppf \\ ϵ_{λ_{2}}^{HbR} ppf & ϵ_{λ_{2}}^{HbO} ppf \end{matrix}) (\begin{matrix} Δ HbR \\ Δ HbO \end{matrix}),

where

ϵ

is the extinction coefficient and ppf is the partial pathlength factor.⁶⁴^,⁶⁵

Then we construct the GLM on the Hb signals:

Eq. (15)

Δ Hb = [\begin{matrix} T & SS & D \end{matrix}] \cdot β + ε,

where

T

is the convolution of modified gamma function and the stimulus sequence as in Eq. (3), SS is the short-separation regressor,

D

is the drift regressor matrix,

β

is the unknown coefficients assigned to each regressor, and

ε

is the residual term. Here

β

are the model coefficients to be estimated and can be approximated as

\hat{β}

by solving the inverse problem:

Eq. (16)

J = [\begin{matrix} T & SS & D \end{matrix}],

Eq. (17)

\hat{β} = {(J^{T} J)}^{- 1} J^{T} Δ Hb .

Then we derive the HRF function:

Eq. (18)

HRF = t \cdot {\hat{β}}^{'},

where

{\hat{β}}^{'}

is a subset of

\hat{β}

that only contains the coefficients assigned to the regressors

T

.

t

is the modified gamma function defined in Eq. (2). We convert the HRF to

Δ OD

using the modified Beer–Lambert law [Eq. (14)]. The

Δ OD

is the input to the image reconstruction step:

Eq. (19)

(\begin{matrix} {Δ OD}_{λ_{1}} \\ {Δ OD}_{λ_{2}} \end{matrix}) = A (\begin{matrix} Δ HbO \\ Δ HbR \end{matrix}),

where

A

is the sensitivity matrix in Eq. (5) that we derived from Monte Carlo simulation (see Sec. 2.3).

Δ HbO

and

Δ HbR

, the two unknowns in this model, are the Hb values at each vertex in the 3D head model. Equation (19) is rewritten as

Eq. (20)

y = A \cdot b,

which is solved as an inverse problem:

Eq. (21)

b = A^{- 1} \cdot y,

where

A^{- 1}

is the Moore–Penrose pseudoinverse with Tikhonov regularization:

Eq. (22)

A^{- 1} = L^{- 1} {({\hat{A}}^{T} \hat{A} + λ_{1} I)}^{- 1} {\hat{A}}^{T},

where

\hat{A}

is the

A

matrix spatially regularized by the term

L

, as defined in Eqs. (11) and (12). The Tikhonov parameter

λ_{1}

adjusts the balance between high spatial frequency noise and image smoothness and was set to 0.01 times the maximum diagonal value of

A^{T} A

.⁶⁰^,⁶³

Three SS conditions were tested: (i) no SS regression; (ii) average SS; and (iii) fake SS. The descriptions of (ii) and (iii) are in Sec. 2.4.

Two image reconstruction models were tested: (i) head: $A = [A_{brain} A_{non - brain}]$ , where $A_{brain}$ contains vertices only from brain cortex and $A_{non - brain}$ has all vertices that are not within the brain cortex, and (ii) brain only: $A = [A_{brain}]$ .

2.7.

Model Performance Evaluation Metrics

The images are created using the values at the time-to-peak. Time-to-peak is always the same time point as we used a predefined modified gamma function to model the HRF curve. We used two metrics to evaluate the performance of each model: the CBR and $t$ -value. Both metrics are calculated at the time point of the highest HRF amplitude. We calculated CBR for the image of each subject as

Eq. (24)

CBR = \frac{mean (V \cdot I_{ROI})}{std (V \cdot I_{background})},

where

V

is the Hb values in one image,

I_{ROI}

is the mask array of the vertices that are within the region of interest (ROI), and

I_{background}

is the mask array of the vertices that are within the background region. After calculating the CBR value for each image on the subject level, we averaged it over all subjects to derive the mean CBR value as the second evaluation matrix.

The $t$ -value is calculated as

Eq. (23)

t - value = \frac{\bar{V}}{σ / \sqrt{n}},

where

\bar{V}

represents the mean Hb value within the ROI/Hb value within the background region averaged across the subjects,

σ

represents the standard deviation of mean Hb value within the ROI/background region across the subjects, and

n

represents the number of subjects. To define the ROI region and the background region, we input a simulated signal that contains only synthetic HRF (no resting data) into the SS-DOT model and derived a noise-free HbO image. In this image, if the Hb value of a vertex was above 90% of the maximum value, it was defined as within the ROI; if the Hb value of a vertex was below 1% of the maximum value, it was defined as within the background region. The ROI region and the background region are shown in Fig. 4(e).

For the experimental dataset, we further calculated a laterality index (LI):⁶⁶^,⁶⁷

Eq. (25)

LI = (C - I) / (I + C),

where

I

represents the number of vertices that have Hb values larger than half of the maximum activation on the ipsilateral hemisphere and

C

refers to the same value on the contralateral hemisphere. LI was calculated for each subject’s reconstructed image as well as for the group average image. The subject level LI was compared using the Kruskal–Wallis test between the models at the significance level of 0.05. The Kruskal–Wallis test was used instead of one-way ANOVA because the LI data are not normally distributed. The post hoc multiple comparison was corrected by the Bonferroni correction method.⁶⁸

2.8.

Computational Load

The simultaneous model is computationally consuming. The runtime of processing the experimental dataset (12 subjects, 3 runs per subject, and 8 mins recordings per run) on a CPU (Intel i9-9900K) is around 40 s for the sequential models, whereas it is about 62 min for the simultaneous models. We developed a GPU algorithm (described in the Supplementary Material) to accelerate it from 62 to 28 min on NVDIA GeForce GTX 1660 Ti.

3. Experiments

3.1.

Simulation Study Design

To test the performance of the model, we simulated datasets by adding synthetic HRF to experimental resting-state data from 11 subjects (see Sec. 3.2). The model was thoroughly examined under various synthetic HRF amplitudes.

To form the synthetic HRF, we simulated a perturbation of HbO and HbR in the image space in the motor region on the brain cortex as a step function in the space with a blob diameter of 10 mm [see Fig. 2(a)]. The amplitude of HbR was half of HbO but negative. In the temporal domain, the brain response time series data are the convolution of the modified gamma function in Eq. (2) [the shape is presented in Fig. 2(b)] and the stimulus sequence (convoluted HRF, abbreviated as cHRF). The stimulus sequences are the same as the sequences used in experimental design but unified as one stimulus condition. Three replications were carried out with different stim marker sequences with randomized ISI. Then we employed the forward matrix $A$ in Eq. (5) to project the perturbation from image space to channel space. In Fig. 2(c), we show an example of the experimental resting state data and the data with cHRF overlayed onto it. We simulated five different HRF amplitudes, from $0.5 \times 10^{- 5}$ to $4.5 \times 10^{- 5} M$ with linear steps of $1 \times 10^{- 5} M$ .

Fig. 2

Simulation of fNIRS data. (a) The perturbation location of both HbO and HbR. The color bar is normalized by the maximum HbO value. (b) The simulated HRF shape. The height of the HRF is represented by HRF_AMP. (c) An example of resting data and resting data overlayed with HRF sequences with three different HRF amplitudes. LPF, low-pass filter of 0.5 Hz.

We tested two different models of SS regressors: (i) average SS: the average OD data over all 19 mm channel time series after first adding the simulated HRF to each channel (see Sec. 3.1) and (ii) fake SS: the average OD data over all 19 mm channels without adding the simulated HRF. The method (ii) is called fake SS because there was no physical 8 mm channel on the probe.

3.2.

Experimental Study Design

Study protocol was approved by the Institutional Review Board of Boston University. All participants provided written informed consent to take part in this study. All participants were healthy young adults ( $n = 12$ ). 11 participants went through one run of 5-min resting state data collection, and all 12 participants went through three 8-min runs of left- and right-hand ball-squeezing task (illustrated in Fig. 3). During the resting state measurement, we asked participants not to move and to remain mentally idle as much as possible for 5 min. During the ball-squeezing task, the participants sat in front of a computer screen (13.3 in. with resolution of $2560 \times 1600$ ) and held one rubber ball (2.5 in. in diameter) in each hand. When the task started, they squeezed the ball in either the left or the right hand according to a visual cue shown on the computer screen accompanied by an audio sound. The stimulus was 5 s, and the rest in between was randomly varied between 5 and 15 s. In each run, there were 15 trials with left-hand stimulus and 15 trials with right-hand stimulus in a fully randomized order. The stimulus presentation was created with PsychoPy software.⁶⁹ fNIRS signals were measured throughout the experiment (two cascaded NIRSport2 $16 \times 16$ devices, NIRx, Berlin, Germany). The probe design was the same layout as in the simulation design [Fig. 1(a)], with two patches covering the motor regions on each hemisphere.

Fig. 3

Experimental paradigm: (a) the 5-min resting state task and (b) the ball-squeezing task.

4. Results

4.1.

Simulation Results

The CBR and $t$ -values at different HRF amplitudes are shown in Fig. 4 for HbO images (see Fig. S2 in the Supplementary Material for HbR) for eight different models: six sequential models, “noSS head,” “avgSS head,” “fakeSS head,” “noSS brain only,” “avgSS brain only,” and “fakeSS brain only,” and two simultaneous models, “avgSS SS-DOT,” “fakeSS SS-DOT.” With increasing simulated HRF amplitudes, CBR shows an increasing trend for all models [Fig. 4(b)]. The fakeSS SS-DOT provides the highest CBR value compared with the other models. However, when the HRF amplitude is small, the improvement is marginal. For sequential models, out of the three SS methods, the “fakeSS” model has the greatest CBR. Out of the two DOT methods, i.e., head and brain-only, whole head reconstruction has higher CBR than brain-only reconstruction. The observation holds true for the HbR CBR values except that noSS head has a higher CBR than fakeSS head with a negligible margin (see Fig. S2B in the Supplementary Material). The error bar is not shown in Fig. 4(b) or Fig. S2B in the Supplementary Material because it would make the visualization difficult. The plots with the error bars can be found in Fig. S9 in the Supplementary Material. The inferential statistics are shown in Tables S7 and S8 in the Supplementary Material.

Fig. 4

Results from the simulated datasets (HbO). (a) An example of the image reconstructed by the fakeSS SS-DOT model (subject #12, HRF $amplitude = 25 μ M$ ). (b) The CBR value under different HRF amplitudes for each model. The average $t$ -value in the (c) ROI region and (d) background region under different HRF amplitudes for each model. (e) The ROI region and the background region.

The $t$ -value in the ROI region [Fig. 4(c)] has an increasing trend with increasing HRF amplitude, whereas the avgSS models decrease at larger HRF amplitudes. The noSS models show the highest $t$ -value except that the fakeSS SS-DOT model exceeds at the largest HRF value. Comparing the two DOT models, the head models always show a higher $t$ -value than brain-only models.

The $t$ -value in background region [Fig. 4(d)] increases for avgSS models, but it is not affected with increasing HRF amplitude for other models. The fakeSS SS-DOT model has the lowest background $t$ -value for HbO but not for HbR (Figs. S2C and S2D in the Supplementary Material).

4.2.

Experimental Results

In Figs. 5 Fig. 6 Fig. 7–8, the group average and individual reconstructed HbO images are presented for the 12 subjects that performed the ball squeezing task. A higher HbO activation is expected on the contralateral hemisphere to the hand performing the ball squeezing task.⁷⁰ In Figs. 5 and 6, we observe such trends in the group level results for all models, except for “No SS brain-only” under the left-hand condition. The LI values for group average images (Table 2) also support that only NoSS brain-only model has a negative LI value.

Fig. 5

Group average HbO images for each model under the condition of left hand. The color bar is normalized by the maximum HbO value.

Fig. 6

Group averaged HbO images for each model under the condition of right hand. The color bar is normalized by the maximum HbO value.

Fig. 7

Subject level HbO images for each model under the condition of left hand. The color bar is normalized by the maximum HbO value.

Fig. 8

Subject level HbO images for each model under the condition of right hand. The color bar is normalized by the maximum HbO value.

Table 2

t-value from subject level images, CBR, and LI value derived from group averaged image for each model on experimental data (HbO). The bold values are the greatest among the models.

HbO		SS-DOT	SS head	noSS head	SS brain-only	noSS brain-only
Left hand	$t$ -value	2.63	2.60	3.21	2.60	3.05
Left hand	CBR	3.10	1.42	1.03	0.55	0.45
	LI	1.00	1.00	1.00	0.71	−0.20
Right hand	$t$ -value	3.00	2.71	3.35	2.64	3.14
Right hand	CBR	6.55	0.87	0.62	0.52	0.39
	LI	1.00	1.00	0.42	0.72	0.14

In the subject-level results, the contralateral hemisphere activation is not obvious for the two brain only models (Figs. 7 and 8). For the noSS head model, the left-hand condition shows such trends but not the right-hand condition. The LI values of the subject-level images (Fig. 9) also indicate the same: the two brain only models have lower LI values (not statistically significant, inferential statistics are shown in Table S3–S6 in the Supplementary Material). Of the four sequential models, the SS head model shows the best results as both the group level and the subject level results show contralateral activation and higher LI values (not statistically significant, inferential statistics are shown in Tables S3–S6 in the Supplementary Material).

Fig. 9

LI of HbO for each model on the experimental dataset for (a) left-hand condition and (b) right-hand condition.

Next, we compare the best sequential model (i.e., SS head model) with the SS-DOT model. The results from SS-DOT are consistent with what we expected both at the group-level and the subject-level. Under the left-hand condition, the laterality is more obvious both on the group level and the subject level compared with “SS head,” but with a larger activation area. The LI values for the SS-DOT model were not significantly different from the SS head model.

The HbR results are also shown in Figs. S3–S6 in the Supplementary Material. The subject-level HbR images were noisier than HbO images, with only group level images showing a meaningful trend. Consistent with the HbO results, the SS head model outperformed the other sequential models. For the SS-DOT model, the group level results also showed lower HbR on the contralateral hemisphere. The LI values were not statistically significantly different across the models (Fig. S7 in the Supplementary Material).

We also calculated CBR and $t$ -values for the experimental images (HbO results are shown in Table 2). The SS-DOT model shows higher CBR values than the sequential models, and the ‘NoSS head’ models show the highest $t$ -values, which is consistent with the simulation results. This can be explained by visually inspecting Figs. 7 and 8, which show that the variance of the images across the subjects is large. A similar trend can also be observed in HbR values, which are presented in Table S2 in the Supplementary Material for the completeness of the results.

The ratio of improvement in CBR versus noSS brain-only for each subject is shown in Fig. 10. With either SS and/or whole head DOT added to the model, the CBR values increase, especially for the subjects who have lower CBR when neither SS nor whole head image reconstruction is implemented. With SS and whole head DOT added to the sequentially model, the CBR is further improved compared with only one of them being implemented. Based on this, the CBR further increases with the SS-DOT model in which SS and whole head DOT are simultaneously implemented.

Fig. 10

Evaluation of CBR for the experimental dataset (12 subjects). All graphs plot the ratio of improvement in SBR versus the model noSS brain-only. HbO, HbR, and HbT are denoted in red circles, blue squares, and green triangles, respectively.

5. Discussion

In this work, we propose a model, SS-DOT, which performs SS regression and DOT image reconstruction methods simultaneously. To test the performance of the new model, we designed a high-density fNIRS probe, and using this probe, we measured fNIRS data from 12 healthy subjects while they were performing a resting-state task and a ball-squeezing task. Based on the resting-state data that we collected, we simulated an fNIRS dataset by overlaying the resting-state data with a synthetic HRF. We also tested two different SS separation models: avgSS uses the average of all channels shorter than 19 mm as the SS regressor after overlaying it with synthetic HRF and fakeSS uses the same averaged signal before the overlay. The avgSS method was motivated by Ref. 8. The fake SS was innovated because we did not have $\sim 8 mm$ SS channels in our probe due to physical limitations. On both simulated and experimental ball-squeezing datasets, we thoroughly tested the SS-DOT model and other sequential models quantitively and qualitatively. From the simulation results, fakeSS SS-DOT showed greater CBR values but only showed marginal benefits at small HRF amplitudes. From the experimental ball-squeezing study, the results are consistent with the simulation study as SS-DOT has the greatest CBR value. The LI values are not significantly different among the models. This indicates the potential of the SS-DOT model to improve the image quality of fNIRS data reconstruction, especially for larger HRF amplitudes. The conventional sequential models are still suitable and practical enough in most cases as they show ability to separate the brain activation from the extracerebral physiological noise in this study.

It was demonstrated that performing SS regression before the DOT image reconstruction could result in an increased contrast-to-noise ratio and decreased image variance across blocks.⁸ We hypothesized that performing SS regression and DOT simultaneously would perform better than performing each separately. Our results show that the SS-DOT model improves the image quality by increasing the CBR value. However, the benefits are marginal at small HRF amplitudes, in which case, the sequential model performs as good as the simultaneous model. This work also presents the performance of other sequential models, thus contributing the information of which SS method and which DOT method would perform better. The no SS model shows the highest $t$ -value, followed by the fakeSS and avgSS models. We speculate that the SS regressor, whether it is 8 mm measurement or average of 19 mm measurements, is prone to removing the brain activation signal, thus lowering the statistical power to detect the brain activation. The whole head DOT models always outperform the brain-only DOT models. Of note, the $t$ -value would be affected by different sample sizes, so one should be careful when comparing our $t$ -values with other studies that have different sample sizes.

Our simulation results showed that CBR goes down for the avgSS models but not for the other models. This is because, when the HRF amplitude increases, the average of the 19 mm SS channels contains more brain activation signals (Fig. S8 in the Supplementary Material). Thus regressing it out results in removing the brain signals. With smaller HRF amplitudes, regressing it also removes part of the brain signals; however, because the brain signal is relatively small compared with large HRF amplitude, it does not affect CBR values as much. Second, we also acquired the experimental fNIRS dataset using ball squeezing task to examine the performance of the sequential and simultaneous models. The results from the sequential models are comparable to Gregg et al.’s work.⁸ However, there are differences in methodologies between the two studies. For example, SS was not regressed out within a GLM framework but rather was subtracted from the long channel signal using the correlation coefficient between the SS and long channel signals as weight. Moreover, HRF was obtained using block averaging. Thus one of the evaluation metrics contrast-to-noise ratio, which was defined as peak response over standard deviation in the prestimulus baseline, could not be employed in this work. This is because our models only result in the modified gamma function, which has no variation in the pre-stimulus baseline. Despite these differences, performing SS regression and then applying DOT has the best performance among the sequential methods, supporting the findings from Gregg et al.⁸

Because we did not have $\sim 8 mm$ SS channels in our probe due to physical limitations, we simulated a “fake” SS channel in our simulation dataset by averaging the resting state data over all of the shorter channels (i.e., 19 mm channels). Because it was only resting state data without overlaying any HRFs, the signal should have not contained any stimulus driven brain activation signals. Thus regressing it out should not present a risk of removing any brain signals, which is supported by the results in Fig. 4(b). However, there is a concern that the non-brain region could also be activated physiologically by the stimulus, which was not presented in this fake channel. Thus this fake channel does not perfectly simulate the real SS data. In the future, with miniaturized optode designs, SS detectors could be implemented in high-density probes, which will allow us to investigate this further.

Future directions include the regularization methods and full-width at half-maximum analysis. The conventional Tikhonov regularization with varying depth sensitivity method was adopted. However, exploratory future research will focus on whether the added information in the SS-DOT model in the matrix before inversion changes singular value spectrum and different regularization parameters or regularization methods would yield better results.

Solving the SS-DOT model, especially the matrix multiplication step, is computationally demanding on a CPU and thus not practical in real-world settings. We took advantage of the GPU to reduce the whole processing time. Despite the GPU acceleration, the job can become even more demanding with more stimulation conditions, number of channels, longer time of measurement, and larger brain coverage. However, we foresee that, in the future, the fast development of GPU hardware and algorithms will resolve this issue.

6. Conclusions

The SS-DOT model improves the image quality by increasing the CBR values by threefold compared with sequential models. However, the improvement is marginal when the HRF amplitude is small. The sequential model still has a strong ability to separate the brain signal from the extracerebral physiological noises if performed with SS and whole head DOT.

Disclosures

Alexander von Lühmann is an employee of NIRx Medical Technologies, a company that sells fNIRS imagers. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as potential conflicts of interest.

Code and Data Availability

The code and data associated with this manuscript will be offered on the bfNIRS cloud service⁷¹ after acceptance.

Acknowledgment

This work was supported by the NIH BRAIN Initiative (Grant No. 1U01EB029856-01).

References

1.

A. Villringer and U. Dirnagl, “Coupling of brain activity and cerebral blood flow: basis of functional neuroimaging,” Cerebrovasc. Brain Metab. Rev., 7 (3), 240 –276 CEMREV 1040-8827 (1995). Google Scholar

2.

F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science, 198 (4323), 1264 –1267 https://doi.org/10.1126/science.929199 (1977). Google Scholar

3.

D. A. Boas et al., “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” Neuroimage, 13 (1), 76 –90 https://doi.org/10.1006/nimg.2000.0674 NEIMEF 1053-8119 (2001). Google Scholar

4.

M. A. Franceschini et al., “Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging,” Psychophysiology, 40 (4), 548 –560 https://doi.org/10.1111/1469-8986.00057 PSPHAF 0048-5772 (2003). Google Scholar

5.

F. Scholkmann et al., “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” Neuroimage, 85 6 –27 https://doi.org/10.1016/j.neuroimage.2013.05.004 NEIMEF 1053-8119 (2014). Google Scholar

6.

D. A. Boas, A. M. Dale and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage, 23 (suppl 1), S275 –S288 https://doi.org/10.1016/j.neuroimage.2004.07.011 NEIMEF 1053-8119 (2004). Google Scholar

7.

D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt., 44 (10), 1957 –1968 https://doi.org/10.1364/AO.44.001957 APOPAI 0003-6935 (2005). Google Scholar

8.

N. Gregg et al., “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front. Neuroenergetics, 2 (July), 1 –8 https://doi.org/10.3389/fnene.2010.00014 (2010). Google Scholar

9.

M. D. Wheelock, J. P. Culver and A. T. Eggebrecht, “High-density diffuse optical tomography for imaging human brain function,” Rev. Sci. Instrum., 90 (5), 1 –24 https://doi.org/10.1063/1.5086809 RSINAK 0034-6748 (2019). Google Scholar

10.

L. Gagnon et al., “Short separation channel location impacts the performance of short channel regression in NIRS,” Neuroimage, 59 (3), 2518 –2528 https://doi.org/10.1016/j.neuroimage.2011.08.095 NEIMEF 1053-8119 (2012). Google Scholar

11.

G. Jasdzewski et al., “Differences in the hemodynamic response to event-related motor and visual paradigms as measured by near-infrared spectroscopy,” Neuroimage, 20 (1), 479 –488 https://doi.org/10.1016/S1053-8119(03)00311-2 NEIMEF 1053-8119 (2003). Google Scholar

12.

M. A. Yücel et al., “Mayer waves reduce the accuracy of estimated hemodynamic response functions in functional near-infrared spectroscopy,” Biomed. Opt. Express, 7 (8), 3078 https://doi.org/10.1364/BOE.7.003078 BOEICL 2156-7085 (2016). Google Scholar

13.

G. Gratton and P. M. Corballis, “Removing the heart from the brain: compensation for the pulse artifact in the photon migration signal,” Psychophysiology, 32 (3), 292 –299 https://doi.org/10.1111/j.1469-8986.1995.tb02958.x PSPHAF 0048-5772 (1995). Google Scholar

14.

G. Morren et al., “Detection of fast neuronal signals in the motor cortex from functional near infrared spectroscopy measurements using independent component analysis,” Med. Biol. Eng. Comput., 42 (1), 92 –99 https://doi.org/10.1007/BF02351016 MBECDY 0140-0118 (2004). Google Scholar

15.

I. Tachtsidis et al., “Functional optical topography analysis using statistical parametric mapping (SPM) methodology with and without physiological confounds,” Oxygen Transport to Tissue XXXI, 237 –243 Springer US( (2010). Google Scholar

16.

S. G. Diamond et al., “Physiological system identification with the Kalman filter in diffuse optical tomography,” Lect. Notes Comput. Sci., 3750 649 –656 https://doi.org/10.1007/11566489_80 LNCSD9 0302-9743 (2005). Google Scholar

17.

S. G. Diamond et al., “Dynamic physiological modeling for functional diffuse optical tomography,” Neuroimage, 30 (1), 88 –101 https://doi.org/10.1016/j.neuroimage.2005.09.016 NEIMEF 1053-8119 (2006). Google Scholar

18.

E. Kirilina et al., “The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy,” Neuroimage, 61 (1), 70 –81 https://doi.org/10.1016/j.neuroimage.2012.02.074 NEIMEF 1053-8119 (2012). Google Scholar

19.

Y. Zhang et al., “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt., 10 (1), 011014 https://doi.org/10.1117/1.1852552 JBOPFO 1083-3668 (2005). Google Scholar

20.

M. A. Franceschini et al., “Diffuse optical imaging of the whole head,” J. Biomed. Opt., 11 (5), 054007 https://doi.org/10.1117/1.2363365 JBOPFO 1083-3668 (2006). Google Scholar

21.

S. Patel et al., “Quantification of systemic interference in optical topography data during frontal lobe and motor cortex activation: an independent component analysis,” Oxygen Transport to Tissue XXXII, 45 –51 Springer US( (2011). Google Scholar

22.

R. Saager and A. Berger, “Measurement of layer-like hemodynamic trends in scalp and cortex: implications for physiological baseline suppression in functional near-infrared spectroscopy,” J. Biomed. Opt., 13 (3), 034017 https://doi.org/10.1117/1.2940587 JBOPFO 1083-3668 (2008). Google Scholar

23.

Q. Zhang, E. N. Brown and G. E. Strangman, “Adaptive filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study,” J. Biomed. Opt., 12 (4), 044014 https://doi.org/10.1117/1.2754714 JBOPFO 1083-3668 (2007). Google Scholar

24.

T. Yamada, S. Umeyama and K. Matsuda, “Multidistance probe arrangement to eliminate artifacts in functional near-infrared spectroscopy,” J. Biomed. Opt., 14 (6), 064034 https://doi.org/10.1117/1.3275469 JBOPFO 1083-3668 (2009). Google Scholar

25.

S. Umeyama and T. Yamada, “Monte Carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy,” J. Biomed. Opt., 14 (6), 064025 https://doi.org/10.1117/1.3275466 JBOPFO 1083-3668 (2009). Google Scholar

26.

N. Gregg et al., “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front. Neuroenergetics, 2 14 https://doi.org/10.3389/fnene.2010.00014 (2010). Google Scholar

27.

L. Gagnon et al., “Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling,” Neuroimage, 56 (3), 1362 –1371 https://doi.org/10.1016/j.neuroimage.2011.03.001 NEIMEF 1053-8119 (2011). Google Scholar

28.

R. B. Saager and A. J. Berger, “Direct characterization and removal of interfering absorption trends in two-layer turbid media,” J. Opt. Soc. Am. A, 22 (9), 1874 https://doi.org/10.1364/JOSAA.22.001874 JOAOD6 0740-3232 (2005). Google Scholar

29.

A. von Lühmann et al., “Improved physiological noise regression in fNIRS: a multimodal extension of the general linear model using temporally embedded canonical correlation analysis,” Neuroimage, 208 116472 https://doi.org/10.1016/j.neuroimage.2019.116472 NEIMEF 1053-8119 (2020). Google Scholar

30.

A. Ortega-Martinez et al., “Multivariate Kalman filter regression of confounding physiological signals for real-time classification of fNIRS data,” Neurophotonics, 9 (2), 025003 https://doi.org/10.1117/1.NPh.9.2.025003 (2022). Google Scholar

31.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 (2), R41 https://doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 (1999). Google Scholar

32.

A. Y. Bluestone et al., “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express, 9 (6), 272 https://doi.org/10.1364/OE.9.000272 OPEXFF 1094-4087 (2001). Google Scholar

33.

A. Custo et al., “Anatomical atlas-guided diffuse optical tomography of brain activation,” Neuroimage, 49 (1), 561 –567 https://doi.org/10.1016/j.neuroimage.2009.07.033 NEIMEF 1053-8119 (2010). Google Scholar

34.

D. K. Joseph et al., “Design and validation of a diffuse optical tomography system with improved spatial resolution for functional brain imaging,” in Opt. InfoBase Conf. Pap., (2006). https://doi.org/10.1364/BIO.2006.WD10 Google Scholar

35.

B. R. White, “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt., 15 (2), 026006 https://doi.org/10.1117/1.3368999 JBOPFO 1083-3668 (2010). Google Scholar

36.

S. R. Arridge et al., “A finite element approach for modeling photon transport in tissue,” Med. Phys., 20 (2), 299 –309 https://doi.org/10.1118/1.597069 MPHYA6 0094-2405 (1993). Google Scholar

37.

E. Okada et al., “Experimental validation of Monte Carlo and finite-element methods for the estimation of the optical path length in inhomogeneous tissue,” Appl. Opt., 35 (19), 3362 https://doi.org/10.1364/AO.35.003362 APOPAI 0003-6935 (1996). Google Scholar

38.

K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys., 22 (6), 691 –701 https://doi.org/10.1118/1.597488 MPHYA6 0094-2405 (1995). Google Scholar

39.

R. Graaff et al., “Condensed Monte Carlo simulations for the description of light transport,” Appl. Opt., 32 (4), 426 –434 https://doi.org/10.1364/AO.32.000426 APOPAI 0003-6935 (1993). Google Scholar

40.

M. Hiraoka et al., “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol., 38 (12), 1859 –1876 https://doi.org/10.1088/0031-9155/38/12/011 PHMBA7 0031-9155 (1993). Google Scholar

41.

S. L. Jacques, L. Wang, “Monte Carlo modeling of light transport in tissues,” Optical Response Laser-Irradiated Tissue, 2607 73 –100 Springer( (1995). Google Scholar

42.

D. A. Boas et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express, 10 (3), 159 –170 https://doi.org/10.1364/OE.10.000159 OPEXFF 1094-4087 (2002). Google Scholar

43.

H. Dehghani et al., “Depth sensitivity analysis of high-density imaging arrays for mapping brain function with diffuse optical tomography,” SPIE Proc., 7174 71741J https://doi.org/10.1117/12.807746 (2009). Google Scholar

44.

Y. Zhan et al., “Singular value decomposition based regularization prior to spectral mixing improves crosstalk in dynamic imaging using spectral diffuse optical tomography,” Biomed. Opt. Express, 3 (9), 2036 https://doi.org/10.1364/BOE.3.002036 BOEICL 2156-7085 (2012). Google Scholar

45.

C. Habermehl et al., “Optimizing the regularization for image reconstruction of cerebral diffuse optical tomography,” J. Biomed. Opt., 19 (9), 096006 https://doi.org/10.1117/1.JBO.19.9.096006 JBOPFO 1083-3668 (2014). Google Scholar

46.

T. J. Huppert et al., “HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Appl. Opt., 48 (10), D280 –D298 https://doi.org/10.1364/AO.48.00D280 APOPAI 0003-6935 (2009). Google Scholar

47.

C. M. Aasted et al., “Anatomical guidance for functional near-infrared spectroscopy: AtlasViewer tutorial,” Neurophotonics, 2 (2), 020801 https://doi.org/10.1117/1.NPh.2.2.020801 (2015). Google Scholar

48.

M. S. Hassanpour et al., “Statistical analysis of high density diffuse optical tomography,” Neuroimage, 85 104 –116 https://doi.org/10.1016/j.neuroimage.2013.05.105 NEIMEF 1053-8119 (2014). Google Scholar

49.

T. J. Huppert, S. G. Diamond and D. A. Boas, “Direct estimation of evoked hemoglobin changes by multimodality fusion imaging,” J. Biomed. Opt., 13 (5), 054031 https://doi.org/10.1117/1.2976432 JBOPFO 1083-3668 (2008). Google Scholar

50.

M. A. Yücel et al., “Calibrating the BOLD signal during a motor task using an extended fusion model incorporating DOT, BOLD and ASL data,” Neuroimage, 61 (4), 1268 –1276 https://doi.org/10.1016/j.neuroimage.2012.04.036 NEIMEF 1053-8119 (2012). Google Scholar

51.

B. B. Zimmermann et al., “Development of a wearable fNIRS system using modular electronic optodes for scalability,” in Biophotonics Congr., BW1A.3 (2019). https://doi.org/10.1364/BRAIN.2019.BW1A.3 Google Scholar

52.

W. D. Penny et al., Statistical Parametric Mapping: The Analysis of Functional Brain Images, Elsevier( (2011). Google Scholar

53.

M. Li et al., “Characterization of the hemodynamic response function in white matter tracts for event-related fMRI,” Nat. Commun., 10 1140 https://doi.org/10.1038/s41467-019-09076-2 NCAOBW 2041-1723 (2019). Google Scholar

54.

Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express, 17 (22), 20178 https://doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 (2009). Google Scholar

55.

D. L. Collins et al., “Design and construction of a realistic digital brain phantom,” IEEE Trans. Med. Imaging, 17 (3), 463 –468 https://doi.org/10.1109/42.712135 ITMID4 0278-0062 (1998). Google Scholar

56.

K. J. Friston, “Statistical parametric mapping,” Neuroscience Databases: A Practical Guide, 237 –250 Springer US( (2003). Google Scholar

57.

G. Strangman, M. A. Franceschini and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage, 18 (4), 865 –879 https://doi.org/10.1016/S1053-8119(03)00021-1 NEIMEF 1053-8119 (2003). Google Scholar

58.

A. Custo et al., “Effective scattering coefficient of the cerebral spinal fluid in adult head models for diffuse optical imaging,” Appl. Opt., 45 (19), 4747 –4755 https://doi.org/10.1364/AO.45.004747 APOPAI 0003-6935 (2006). Google Scholar

59.

S. L. Jacques, “Erratum: optical properties of biological tissues: a review,” Phys. Med. Biol., 58 (14), 5007 –5008 https://doi.org/10.1088/0031-9155/58/14/5007 PHMBA7 0031-9155 (2013). Google Scholar

60.

W. Fan, H. Dehghani and A. T. Eggebrecht, “Investigation of effect of modulation frequency on high-density diffuse optical tomography image quality,” in Opt. InfoBase Conf., 1 –19 (2021). https://doi.org/10.1364/ECBO.2021.EM3C.4 Google Scholar

61.

S. Wray et al., “Characterization of the near infrared absorption spectra of cytochrome aa3 and haemoglobin for the non-invasive monitoring of cerebral oxygenation,” Biochim. Biophys. Acta, 933 (1), 184 –192 https://doi.org/10.1016/0005-2728(88)90069-2 (1988). Google Scholar

62.

S. Jahani et al., “Motion artifact detection and correction in functional near-infrared spectroscopy: a new hybrid method based on spline interpolation method and Savitzky–Golay filtering,” Neurophotonics, 5 (1), 015003 https://doi.org/10.1117/1.NPh.5.1.015003 (2018). Google Scholar

63.

J. P. Culver et al., “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys., 30 (2), 235 –247 https://doi.org/10.1118/1.1534109 MPHYA6 0094-2405 (2003). Google Scholar

64.

D. T. Delpy et al., “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol., 33 (12), 1433 –1442 https://doi.org/10.1088/0031-9155/33/12/008 PHMBA7 0031-9155 (1988). Google Scholar

65.

S. J. Matcher et al., “Performance comparison of several published tissue near-infrared spectroscopy algorithms,” Anal. Biochem., 227 (1), 54 –68 https://doi.org/10.1006/abio.1995.1252 ANBCA2 0003-2697 (1995). Google Scholar

66.

S. N. Maxfield Panker, “The effects of robotic training and cortical stimulation on reaching skill after chronic stroke,” (2012). Google Scholar

67.

Q. Tang et al., “Modulation of interhemispheric activation balance in motor-related areas of stroke patients with motor recovery: systematic review and meta-analysis of fMRI studies,” Neurosci. Biobehav. Rev., 57 392 –400 https://doi.org/10.1016/j.neubiorev.2015.09.003 NBREDE 0149-7634 (2015). Google Scholar

68.

O. J. Dunn, “Multiple comparisons among means,” J. Am. Stat. Assoc., 56 (293), 52 –64 https://doi.org/10.1080/01621459.1961.10482090 (1961). Google Scholar

69.

J. W. Peirce, “Generating stimuli for neuroscience using PsychoPy,” Front. Neuroinform., 2 (JAN), 1 –8 https://doi.org/10.3389/neuro.11.010.2008 (2009). Google Scholar

70.

P. Lapborisuth et al., “Neurofeedback-based functional near-infrared spectroscopy upregulates motor cortex activity in imagined motor tasks,” Neurophotonics, 4 (2), 021107 https://doi.org/10.1117/1.NPh.4.2.021107 (2017). Google Scholar

71.

Y. Gao, “BallSqueezingHD,” https://bfnirs.openfnirs.org/ (). Google Scholar

Biography

Yuanyuan Gao received her bachelor’s and master’s degrees from Beihang University and her PhD from Rensselaer Polytechnic Institute. She studied the effects of neuromodulation on human motor behaviors using fNIRS for her PhD research. She worked as a postdoctoral fellow at Boston University Neurophotonics Center. Her research interests include functional near-infrared spectroscopy (fNIRS), multimodeled with electroencephalography (EEG), eye-tracking, physiology measurements, neuromodulation, and machine-learning models. She started her second term of postdoctoral training at Stanford University starting in March 2023.

De’Ja Rogers received her BS degree in electrical and electronics engineering from Norfolk State University and her MS degree in biomedical engineering from Boston University. She is a biomedical engineering PhD student at the Bio Optical Acoustic Spectroscopy Lab of Boston University. She uses the combination fNIRS and EEG to measure hemodynamic oscillations during all sleep stages, with the goal of investigating sleep’s relation to neurodegeneration and validating fNIRS as a sleep study instrument.

Alexander von Lühmann received his BSc and MSc degrees in electrical engineering from Karlsruhe Institute of Technology in 2011 and 2014, respectively, and his PhD (Dr-Ing) with distinction in 2018 from Technische Universität (TU) Berlin. He is a head of the “Intelligent Biomedical Sensing” group in the Machine Learning Department of TU Berlin–Berlin Institute for the Foundations of Learning and Data. Formerly, he was an R&D director at NIRx Medizintechnik GmbH, a postdoc at Boston University, and a visiting researcher at Harvard Medical School.

Antonio Ortega-Martinez received his PhD from the Bio Optical Acoustic Spectroscopy Lab at Boston University and an electrical engineer with a master’s degree in optics. His career has most often revolved around industrial and biomedical applications of optics and infrared technologies. He performed research on wearable fNIRS devices, signal processing, brain–computer interfaces, and time-domain fNIRS.

David A. Boas is a professor of biomedical engineering at Boston University. He was the founding president of the Society for Functional Near-Infrared Spectroscopy and founding Editor-in-Chief of the SPIE journal Neurophotonics. He received the Britton Chance Biomedical Optics Award in 2016 for his development of several innovative high-impact biomedical optical technologies in the neurosciences, as well as following through with impactful application studies, and fostering the widespread adoption of these technologies.

Meryem Ayşe Yücel is a research associate professor at Boston University (BU). Prior to her position at BU, she was an assistant in biomedical engineering at Massachusetts General Hospital and an instructor at Harvard Medical School, Radiology. Her primary research interest is to understand how the brain works in health and disease. Throughout her career, she has gained expertise in mathematical modeling of biological systems and functional brain imaging (fNIRS, fMRI, and EEG).

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Yuanyuan Gao, De’Ja Rogers, Alexander von Lühmann, Antonio Ortega-Martinez, David A. Boas, and Meryem A. Yücel "Short-separation regression incorporated diffuse optical tomography image reconstruction modeling for high-density functional near-infrared spectroscopy," Neurophotonics 10(2), 025007 (23 May 2023). https://doi.org/10.1117/1.NPh.10.2.025007

Received: 14 November 2022; Accepted: 3 May 2023; Published: 23 May 2023

Access the abstract

JOURNAL ARTICLE
18 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY