Open Access
Add to BibSonomy
Abstract
High density diffuse optical tomography has become increasingly important to detect underlying neuronal activities. Conventional methods first estimate the time courses of the changes in the absorption coefficients for all the voxels, and then estimate the hemodynamic response function (HRF). Activation-level maps are extracted at last based on this HRF. However, the error propagation among the successive processes degrades and even misleads the final results. Besides, the computation burden is heavy. To address the above problems, a direct method is proposed in this paper to simultaneously estimate the HRF and the activation-level maps from the boundary fluxes. It is assumed that all the voxels in the same activated brain region share the same HRF but differ in the activation levels, and no prior information is imposed on the specific shape of the HRF. The dynamic simulation and phantom experiments demonstrate that the proposed method outperforms the conventional one in terms of the estimation accuracy and computation speed.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Optical neuroimaging techniques have shown promise to understand the brain function because of the advantages including good temporal resolution, portability, quietness and low sensitivity to motion artifacts [1,2]. In recent years, they have evolved into high density diffuse optical tomography (HD-DOT) from the initial functional near-infrared spectroscopy. They both can distinguish the concentration changes in oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR), which are directly associated with the metabolisms of activated brain regions [3–5]. However, different from the functional near-infrared spectroscopy, HD-DOT utilizes the overlapping measurements at multiple distances from dense regular arrays of sources and detectors. Besides, HD-DOT inverts the physically consistent photon-migration-model to reconstruct a 3-dimensional image and significantly improves the spatial resolution and quantitative accuracy.
In the conventional HD-DOT, the changes in absorption coefficients (${\mu _a}$) at every time point are separately reconstructed and a time course for every voxel is formed. Then the changes in ${\mu _a}$ before and after the stimulation onset are used to help judge whether a brain region is activated. However, the HD-DOT data is usually contaminated by the physiological interferences induced by cardiac pulsations, respiratory and blood pressure changes [6]. These physiological interferences may distort the time courses and even mislead the judgment of the activated regions. To address this problem, the general linear model (GLM) is introduced to explain the time course of the ${\mu _a}$ changes as the linear combination of some explanatory variables [7,8]. The coefficients of the convolution of stimulus onsets with the hemodynamic response function (HRF) are referred to as the activation levels. They are more reliable to determine whether a brain region is activated.
In the above GLM method, the HRF accuracy is important for the estimation of the activation levels. Like functional magnetic resonance imaging (fMRI), in DOT the brain is also treated as a linear time invariant system [9]. The HRF is defined as the impulse response evoked by a short and single stimulus of unit intensity [10]. The power of statistical inferences based on the GLM will also be weakened by the inaccurate HRF. Furthermore, estimating the precise timing and waveform of the HRF is meaningful to explore the relative timing of neuronal activities, neuronal feedback processes and sustained activities within a brain region [11].
For estimating the HRF, there have been many methods proposed in fMRI and DOT [12]. These methods vary greatly in the degree to which they make a priori assumptions about the shape of the HRF [13]. At one extreme, the shape of the HRF is assumed to be pre-determined [13,14]. Then, some methods are proposed to assume the HRF to be the Poisson, Gamma, Gaussian and joint half cosine functions. However, experiments have shown that the HRF varies across different subjects, ages, days, stimuli and brain regions [15]. Thus, it is necessary to devise more feasible schemes to estimate the various shapes of HRF. Then the HRF is assumed to be a linear combination of some basis function sets. The most common choice is to use the canonical HRF and its derivatives with respect to time and dispersion as the basis function sets [14,16]. In addition, some methods based on cosine functions, the principal components, radial basis functions, spline basis sets and spectral basis function are proposed [13]. Furthermore, the most flexible method assumes no a priori knowledge about the HRF shape. For example, the finite impulse response model has been proposed within the GLM framework [17]. It contains one free parameter for every time point following the stimulation onsets. However, this flexibility has costs. The degree of freedom is too large and it is sensitive to the measurement noise. Thus, the temporal regularization has been introduced in the smooth finite impulse response (SFIR) to favor the HRF with small second order temporal derivatives [18].
However, all the above methods have two limitations. First, the HRF and the activation levels are indirectly estimated from the time courses of ${\mu _a}$ changes instead of being directly estimated from the boundary fluxes. The error in the separate computation of the ${\mu _a}$ changes at all the time points will transfer to the subsequent estimations of the activation-level maps. If the estimated HRF is not accurate, the activation-level maps will be degraded and even misled. Furthermore, both the HRF and the activation levels are determined based on the separate time courses of each voxel instead of combining the time course of all the voxels. When the signal to noise ratio (SNR) is low, the accurate estimation is difficult. To overcome these shortcomings, a joint direct estimation (JDE) method is proposed to simultaneously estimate the HRF and the activation-level maps from the boundary fluxes. Based on the low-rank constraint and the assumption that all the voxels in the same activated region share the same HRF, it avoids solving the ${\mu _a}$ changes for every time point. Besides, it imposes no constraint on the specific shape of the HRF. Furthermore, the a priori smooth information is utilized to improve the estimated HRF.
2. Methods
The inverse problems in DOT at all the time points (k=1, 2, ……, K) can be stacked together as follows:
where ${\bf M} = [{{{\bf m}_1}{{\bf m}_2} \cdots {{\bf m}_K}} ]\in {R^{M \times K}}$ with ${{\bf m}_k}$ being the boundary flux changes at the k-th time point and $M$ being the total number of the measurements; ${\bf X} = [{{{\bf x}_1}{{\bf x}_2} \cdots {{\bf x}_K}} ]\in {R^{N \times K}}$ with ${{\bf x}_k}$ being the ${\mu _a}$ changes at the k-th time point and N being the total number of the voxels; ${\bf J} \in {R^{M \times N}}$ is the Jacobian matrix to describe the sensitivity of ${{\bf m}_k}$ to ${{\bf x}_k}$.According to the GLM, ${\bf X}$ can be expressed in the following form:
Two independent random variables are orthogonal when at least one of the random variables has zero mean (see Appendix A) [19]. The nuisance interference signals in the scalp (${\bf f}$) and hemodynamic signals in the brain (${\bf h}$) are usually assumed to be independent from each other [20,21]. ${\bf f}$ denotes the periodic changes of ${\mu _a}$ in the scalp caused by the cardiac, respiration and Mayer wave. It can be usually simulated by sinusoidal functions [20]. Thus, its mean value is 0. ${\bf h}$ and ${\bf f}$ can be further inferred to be orthogonal each other.
Suppose that $r({\cdot} )$ denote the rank of a vector or matrix. $r({\boldsymbol{\mathrm{\beta}} })$, $r({\bf h})$, $r({\boldsymbol{\mathrm{\alpha}} })$ and $r({\bf f})$ are all equal to 1. The product of multiple matrices is limited in its rank to the lowest of the constituent matrix ranks (see Appendix B) [22]. Thus $r({\bf {J}\boldsymbol{\mathrm{\beta}} }{({\bf {Sh}})^T}) = 1$ and $r({\bf {J}\boldsymbol{\mathrm{\alpha}} }{{\bf f}^T}) = 1$. On the one hand, the rank of the sum of two matrices is smaller than or equal to the sum of the ranks of the two matrices, thus $r({\bf {M}}) \le 2$. On the other hand, because ${\bf M}$ has at least two independent components ${\bf h}$ and ${\bf f}$, $r({\bf M}) \ge 2$ [23,24]. Therefore, we obtain $r({\bf M}) = 2$.
In the singular value decomposition of ${\bf M}$, the singular vectors are orthogonal each other and the singular values are sorted in a descending order. The bigger singular value, the bigger weight of the singular vector. Since $r({\bf M}) = 2$, the problem in Eq. (3) is essentially to find a rank-2 approximation to the measurement ${\bf M}$. According to the Eckart and Young theorem, a rank-2 approximation can be obtained by the truncated singular value decomposition method [25]. Specifically, ${\bf M}$ can be approximated by
where ${\xi _i}$, ${{\bf u}_i}$ and ${\bf v}_i^{}$ are respectively the i-th singular value and left singular vector and right singular vector. In addition, the boundary measurements exponentially decrease with increased distance from the scalp. ${\bf M}$ is more sensitive to the scalp interference signals ${\bf f}$ than the cerebral signals ${\bf h}$ induced by the brain activation [26]. Thus, ${\bf f}$ has bigger weight than ${\bf h}$. Therefore, the first term and second term in Eq. (4) respectively account for ${\bf {J}\boldsymbol{\mathrm{\alpha}} }{{\bf f}^T}$ and ${\bf {J}\boldsymbol{\mathrm{\beta}} }{({\bf {Sh}})^T}$. Accordingly, the second term can be expressed as follows: By further decomposing Eq. (5), it is easy to obtain ${\bf {J}\boldsymbol{\mathrm{\beta}} } = {{\bf u}_2}$ and ${\bf {Sh}} = {\xi _2}{{\bf v}_2}$. By solving separately them, ${\boldsymbol{\mathrm{\beta}} }$ and ${\bf h}$ can be solved.In practice, the changes in the HbO and HbR concentrations are induced by the blood volume and flow. Therefore, the true ${\bf h}$ is a smooth curve. This a prior information is helpful to regularize ${\bf h}$. The equation ${\bf{Sh}} = {\xi _2}{{\bf v}_2}$ is converted into the following optimization problem by the commonly used Tikhonov regularization method [22]:
3. Simulation experiments
To investigate the performances of the proposed method, some simulation experiments with block paradigms are conducted. Another commonly used non-parametric method based on the SFIR is used as a reference method to compare with JDE. To quantitatively compare these two methods, the estimated HRF and the true HRF are first normalized as ${{\bf h}_{est}}$ and ${{\bf h}_{tr}}$, respectively. Then, the root mean square error (RMSE) is defined below [22]
3.1 Optical model
When photons are emitted into the brain tissues, they are often absorbed and scattered multiple times before arriving at detectors. The Monte Carlo and the diffusion equation are commonly used to model the transport of the photons in the turbid medium. Considering the low-scattering property of the cerebral spinal fluid (CSF) where the diffusion equation is not applicable, the golden standard method Monte Carlo is used. To accelerate the computation speed, the Monte Carlo based on the graphics processing unit is used to perform the parallel computation [28]. A 5-layered MRI head atlas is utilized and it includes the scalp, skull, CSF, gray matter and white matter [29]. The detailed layered structure of the atlas is shown in Fig. 1. The concentrations of HbO and HbR for all the layers are listed in Table 1 and their optical properties at the wavelengths of 760 nm and 830 nm are detailed in Table 2 [30]. 20 light sources (red circles) and 12 detectors (blue squares) are assigned on the scalp surface with an equal interval of 20 mm as shown in Fig. 2. The coordinates of Y and Z for the sources range from 50 mm to 110 mm and 48 mm to 128 mm, respectively. Because the sources are placed on the scalp, their X coordinates are determined by the outermost scalp with the same Y and Z coordinates with the sources. The detectors are located at the centers of every 4 joint sources. Likewise, their X coordinates are determined by the outermost scalp with the same Y and Z coordinates with the detectors. One source-detector pair is defined as a channel and all the channels are divided into the first, the second and the third nearest-neighbor (NN) channels (NN-1, NN-2 and NN-3). Their source-detector separations are 13.44 mm, 30.04 mm and 40.31 mm, respectively. Because the measurements of NN-2 and NN-3 have appropriate detection depth and SNR, they are used to estimate the HRF and the activation-level maps. Besides, in the following simulation experiments, the voxels in the gray matter in the domain of $24mm < X < 50\textrm{ mm}$, $\textrm{77 mm < Y < 85 mm}$ and $\textrm{80 mm < Z < 100 mm}$ consist the activated brain region, as shown in Fig. 1 and Fig. 2.
Fig. 1. Layered head structure: (a) vertical plane at Z=95 mm; (b) coronal plane at X=40 mm. Different colors denote different tissue types.
Fig. 2. Brain atlas: (a) right view; (b) back view with assignment of NN-1, NN-2 and NN-3. The green color denotes the activated brain region.
Table 1. Concentrations of HbO and HbR for different tissue types.
Table 2. Optical properties for different tissue types.
The HRF of HbO concentration is modeled as the linear combination of two gamma functions [31,32]:
Fig. 3. Normalized time courses: (a) ${h_{HbO}}(k)$ and ${h_{HbR}}(k)$; (b) The stimulus and its convolution with ${h_{HbO}}(k)$ and ${h_{HbR}}(k)$; (c) Changes of ${\mu _a}$ in the activated brain region at the wavelengths of 760 nm and 830 nm; (d) The scalp interference signal.
As shown in Fig. 3(b), the block stimuli are consisted of two box functions whose width and interval are 5 s and 20 s, respectively. The normalized convolution of HRF (HbO and HbR concentrations) and stimulation are also presented in Fig. 3(b). According to the concentrations of HbO and HbR in Table 1, the normalized time course of ${\mu _a}$ changes at the wavelength of 760 nm and 830 nm are calculated and shown in Fig. 3(c). Besides to the changes in ${\mu _a}$ happening in the gray matter, the interference often occurs in the scalp and hinders the estimation of the HRF and the activation-level map. They are usually caused by the cardiac, respiration and Mayer wave, which can be simulated with sinusoidal oscillations with frequencies being 1 Hz, 0.3 Hz and 0.1 Hz, respectively. The normalized waveform of the sum of scalp interferences is shown in Fig. 3(d). The ${\mu _a}$ in the gray matter and scalp both change at most about 50% relative to their own initial values. The total temporal length is 60 s and the temporal sampling frequency is 10 Hz. Therefore, in total 600 time points are considered.
Though some researchers assume the scalp interferences to be global, it is still necessary to investigate the local scalp interference in the following simulation experiments [34–37]. For the local case, the interferences are assumed to happen only in the scalp in the domain $24mm < X < 50\textrm{ mm}$, $\textrm{50 mm < Y < 100 mm}$ and $\textrm{40 mm < Z < 75 mm}$.
The NN-1 channels with a short source-detector separation are only sensitive to the ${\mu _a}$ changes in the scalp. They can be used as reference signals in the SFIR to filter out the interference components included in the time course of the ${\mu _a}$ changes in the gray matters [38]. Furthermore, the frequency of the cardiac is about 1 Hz, not overlapped by that of the HRF. Thus, a Butterworth-type low pass filter with cut-off frequency of 0.4 Hz is used to remove the high frequency components. At last, a moving average filter with a window width of 5 is used to further smooth the time courses. Besides, the SNR of the measurements is proportional to the squared root of the light intensity at each channel [39]. The Gaussian white noise is added into the pure simulated data to achieve SNR=20 dB for the channels with the weakest light intensity. After the simulation experiments are repeated 10 times, the average values and the standard deviations are used to assess the performances of the proposed method.
3.2 Global scalp interference
Taking the wavelength of 760 nm as an example, the curve of singular values in a decreasing order is shown in Fig. 4(a). It decreases rapidly and the first two singular values are much bigger than the other ones. Besides, the first right singular vector shown in Fig. 4(b) is consistent with the true scalp interference. Likewise, the second right singular vector in Fig. 4(c) is also in accordance with the true changes of ${\mu _a}$ in the activated brain region. These observations confirm the conclusions in Section 2.
Fig. 4. Normalized time courses in the case of global scalp interference: (a) singular value curve; the (b) first and (c) second right singular vector of M.
The time courses of changes in ${\mu _a}$ at the two wavelengths are converted into the changes of HbO and HbR concentrations. The HRFs of HbO and HbR concentrations are estimated by SFIR and JDE. As shown in Fig. 5, HRFs of HbO and HbR concentrations of JDE are closer to the true ones, compared to the SFIR. For the HbR of SFIR, the lowest value appears earlier than the true one. JDE obtains smaller RMSEs than the SFIR, as shown in Table 3.
Table 3. RMSEs of the normalized HRFs of HbO and HbR concentrations estimated by SFIR and JDE in the case of global scalp interference
As shown in Fig. 6, the normalized activated regions estimated by JDE have less artifacts, compared with those of SFIR. Furthermore, the metrics of the activation-level maps are calculated and shown in Table 4. Except for the AR at the wavelength of 760 nm, the JDE method obtains the smaller $RMS{E_\beta }$, bigger CNR, AR closer to 1, demonstrating that the JDE has better expressions to estimate the activation-level maps. In addition, the durations of JDE are much shorter than those of the SFIR. The standard deviations of all the metrics of the JDE are smaller than those of the SFIR. This demonstrates that the JDE is more robust to the noise than the SFIR.
Fig. 6. Coronal planes at X=40 mm of the normalized activation-level maps in the case of global interference: (a) 760 nm, SFIR and (b) 760 nm, JDE; (c) 830 nm, SFIR and (d) 830 nm, JDE; (e) True.
Table 4. Metrics of the normalized activation-level maps estimated by SFIR and JDE in the case of global scalp interference
3.3 Local scalp interference
In the case of local scalp interference, we also take the wavelength of 760 nm as an example. The normalized singular value curve is shown in Fig. 7(a) and it decreases rapidly. But the relative value of the second biggest one is bigger than that of the global interference. This is possibly because the local scalp interference contributes less to the boundary measurements than the global interference. Besides, the first component (Fig. 7(b)) and the second component (Fig. 7(c)) respectively account for the scalp interference and the changes of ${\mu _a}$ in gray matter. Then the changes of ${\mu _a}$ are converted into the changes of HbO and HbR concentrations which are used to estimate the HRFs of HbO and HbR concentrations. As presented in Fig. 8, JDE estimates HRFs closer to the true ones. The HRF of HbR concentration has an artificial positive peak. In Table 5, the JDE method also has smaller RMSEs than the SFIR.
Fig. 7. Normalized time courses in the case of local scalp interference: (a) singular value curve; the (b) first and (c) second right singular vector of M.
Table 5. RMSEs of the normalized HRFs of HbO and HbR concentrations estimated by SFIR and JDE in the case of local scalp interference
In Fig. 9, the normalized activation level maps estimated by the JDE have less artifacts than those of SFIR. The metrics of the activation level maps are compared in Table 6. It is clear that JDE almost improves all the metrics, except the AR. In addition, the JDE has much smaller durations than the SFIR.
Fig. 9. Coronal planes at X=40 mm of the normalized activation-level maps in the case of local interference: (a) 760 nm, SFIR and (b) 760 nm, JDE; (c) 830 nm, SFIR and (d) 830 nm, JDE; (e) True.
Table 6. Metrics of the normalized activation-level maps estimated by SFIR and JDE in the case of local scalp interference
4. Phantom experiments
The performances of the proposed method are further investigated based on dynamic phantom experiments [27].
4.1 Experimental setup
As shown in Fig. 10, a cuboid-shaped daicel phantom of size $100\textrm{ mm } \times \textrm{ 120 mm } \times 3\textrm{0 mm}$ is used to simulate a human head. At the wavelength of 785 nm, its optical properties are about ${\mu _a} = 0.0041/mm$, ${\mu ^{\prime}_s} = 1/mm$[40]. Inside this phantom, a cylinder hole is embedded to simulate the activated brain region. Its diameter and depth are 15 mm and 18 mm, respectively. It is worth noting that, the distance between the bottom of the cylinder hole and the bottom of the daicel phantom is 12 mm, which is thick enough to simulate the scalp. The mixed solution of intralipid (mimicking scattering) and ink (mimicking absorption) are used to simulate the brain at different states. Cup A and B are respectively filled with target solution (${\mu _a} = 0.0260/mm$, ${\mu ^{\prime}_s} = 1.0134/mm$) and background solution (${\mu _a} = 0.0173/mm$, ${\mu ^{\prime}_s} = 1.0134/mm$). The mixed solution in Cup A and B is pumped into the cylinder hole by an input constant flow pump. At the same time, the solution in the cylinder hole is pumped into Cup C by another output constant flow pump. Two pumps’ flows are adjusted to be equal so that the volume of the solution in the cylinder hole keeps constant. At the beginning, the cylinder hole is filled with the background solution. After operating the equipment according to Table 7, ${\mu _a}$ in the cylinder hole will first increase and then decrease back to the initial value.
Fig. 10. Schematic of the dynamic phantoms and the HD-DOT detecting equipment. Cup A, B and C are filled with the target solution, background solution, and waste solution, respectively. Switch 1 and 2 are used to select a cup from which the solution is pumped into the cylinder hole.
Table 7. Operation process of the dynamic phantom equipment.
As shown in Fig. 10, a HD-DOT prototype instrument based on a lock-in photon-counting technique is utilized [41]. It combines the high sensitivity of the photon-counting technique and the parallel features of the lock-in technique. 20 laser diode sources and 12 photon counting photomultipliers (H10682, Hamamatsu, Japan) are placed under the phantom. The sources are modulated by square waves with frequencies ranging from 6 kHz to 10 kHz with an interval of 0.2 kHz. The sources are working in a multi-field illumination mode. In every filed, multiple sources are lit simultaneously. The detected photons are demodulated according to the modulating frequencies. In total, 8 fields are needed to light all the sources. The accumulation period for each field is 200 ms and the frame rate is 0.625 Hz. During the whole experimental process, the sources are lit field by field and the boundary fluxes are continuously detected.
4.2 Results
In this phantom experiment, it is very difficult or nearly impossible to accurately know the true HRF of the system. The next best choice is to find a pair of stimulation and HRF whose convolution optimally fit the reconstructed time course of ${\mu _a}$ changes ($\delta {\mu _a}$) in the cylinder hole, as shown in Fig. 11(a). This optimal stimulation is then used in the SFIR and JDE to estimate the HRF, which is then compared with the optimal HRF. In Fig. 11(b), the normalized HRF estimated by JDE is closer to the optimal one than that of SFIR. $RMS{E_h}$ of JDE is smaller than that of SFIR, as shown in the first row in Table 8.
Fig. 11. Normalized time courses of (a) stimulus, HRF and convolution results and (b) HRF estimated by SFIR and JDE.
Table 8. Metrics of the normalized HRF and activation-level maps estimated by SFIR and JDE in the dynamic phantom experiments.
The normalized activation-level maps and their horizonal profiles passing through the peak activation level are plotted in Fig. 12. It is easy to see that the activated region estimated by SFIR is much bigger than the true one. However, the activated region estimated by the JDE is very close to the true size and has less artifacts. The profile of the JDE is narrower than that of SFIR. The metrics of the normalized activation-level maps are compared in Table 8 where the JDE obtains smaller $RMS{E_\beta }$, bigger CNR and AR closer to 1. The last row of Table 8 shows that the JDE has a shorter duration than SFIR.
Fig. 12. Results of the phantom experiments: normalized activation-level maps estimated by (a) SFIR and (b) JDE; The black dotted circles denote the true locations of the activate regions; (c) Horizonal profiles passing through the peak pixel.
5 Discussions
Limited by the existing phantom experimental conditions, the physiological interferences cannot be simulated. Accordingly, the terms associated with the interferences in Eqs. (2), (3) and (4) should be removed. The singular vectors of ${\bf M}$ only include the cerebral signals induced by the brain activation. Therefore, the interference term in Eq. (4) can be neglected and it can be approximated only with the first term:
The activation-level maps and HRF are respectively estimated based on the ${{\bf u}_1}$ and ${\bf v}_1^{}$. Though the phantom experiments cannot completely mimic the real activation processes happening in the human brains, it still to some extent demonstrates the expected better performances of the JDE than the SFIR. Besides, in the phantom experiments, the semi-three-dimensional framework is used to reduce computational burden and improve image qualities. It provides a trade-off between the full 3 dimensional optical tomography and 2 dimensional topography [42].To fairly compare the JDE and SFIR, the stable and fast algebraic reconstruction technique with no a priori information is used when solving the inverse problems in the above sections. However, it is possible to introduce the sparsity regularization to further improve the image qualities [43]. When activated by some stimulus, only some localized region is activated, but the remaining parts are still at rest state. This means that the activation-level map is itself sparse [44]. Whereas, the selection of the optimal regularization parameter is difficult [45]. The commonly used cross validation needs to solve the inverse problems for every optional regularization parameter. This is time consuming especially for the conventional SFIR because it needs to solve one inverse problem to obtain the ${\mu _a}$ changes for every time point. The computation burden will increase with the number of time points. It will be nearly impossible to implement the cross validation for the practically big number of time points. However, the JDE only needs to solve one inverse problem to estimate the activation-level maps. Thus the JDE is much faster than the SFIR, which makes it much easier to introduce the sparsity regularization and select the best regularization parameter. Herein, the well-known L1-norm based fast iterative shrinkage-thresholding algorithm and the regularization parameter of 0.5 are selected to show the improved performances contributed by the sparsity regularization [46]. The activation-map and the profile passing through the peak voxel are presented in Fig. 13. Compared with the results in Fig. 12, the activated region is closer to the true one and the profile is narrower. Besides, the metrics are also improved compared with those in Table 8. $RMS{E_\beta }$, CNR and AR are respectively 0.52, 17.90 and 0.91.
Fig. 13. Results of the phantom experiments: (a) Normalized activation-level maps estimated by the sparsity regularized JDE; (b) Normalized profiles along the horizonal line passing through the peak pixel.
6 Conclusions
By incorporating the low-rank constraint, we propose a direct method to simultaneously estimate the normalized HRF and the activation-level maps from the boundary fluxes. It is unnecessary to separately obtain the ${\mu _a}$ changes at every time point. Besides, the a priori smooth information is imposed on the HRFs of HbO and HbR concentrations to improve the qualities. The better expression of the proposed JDE method is then validated by the dynamic numerical and phantom experiments. The experimental results demonstrate that the JDE can obtain more accurate HRFs. Besides, the activation-level maps estimated by the JDE method have better qualities compared with the conventional SFIR. In conclusion, the proposed JDE method can improve the estimated HRF, activation-level maps and the computation speed.
Appendix A
If random variables A and B are independent, it implies ${p_{X,Y}}(x,y) = {p_X}(x){p_Y}(y)$. Then
If at least one of ${{\cal E}}(X)$ and ${{\cal E}}(Y)$ is equal to 0, ${{\cal E}}(XY) = 0$. According to the definition, random variables X and Y are orthogonal if ${{\cal E}}(XY) = 0$[19].'
Appendix B
The range of an m-by-n matrix A is defined below [47]
Since $r({\bf {AB}}) = r({{\bf B}^T}{{\bf A}^T})$, similar to the above derivation, $r({{\bf B}^T}{{\bf A}^T}) \le r({{\bf B}^T})$[48]. Because $r({{\bf B}^T}) = r({\bf B})$, $r({\bf {AB}})) \le r({\bf B})$. In summary, $r({\bf {AB}})) \le \min (r({\bf A}),r({\bf B})))$ where $\min ({\cdot} )$ denotes the minimization operation.
Funding
National Natural Science Foundation of China (61575140, 81871393, 81971662); Tianjin Municipal Government of China (17JCQNJC12700).
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. S. Lloyd-Fox, A. Blasi, and C. E. Elwell, “Illuminating the developing brain: the past, present and future of functional near infrared spectroscopy,” Neurosci. Biobehav. Rev. 34(3), 269–284 (2010). [CrossRef]
2. L. Gagnon, K. Perdue, D. N. Greve, D. Goldenholz, G. Kaskhedikar, and D. A. Boas, “Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling,” NeuroImage 56(3), 1362–1371 (2011). [CrossRef]
3. A. Shah and A. K. Seghouane, “An integrated framework for joint HRF and drift estimation and HbO/HbR signal improvement in fNIRS data,” IEEE Trans. Med. Imaging 33(11), 2086–2097 (2014). [CrossRef]
4. S. G. Diamond, T. J. Huppert, V. Kolehmainen, M. A. Franceschini, J. P. Kaipio, S. R. Arridge, and D. A. Boas, “Dynamic physiological modeling for functional diffuse optical tomography,” NeuroImage 30(1), 88–101 (2006). [CrossRef]
5. 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), 051101 (2019). [CrossRef]
6. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. Mata Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85(1), 6–27 (2014). [CrossRef]
7. J. Cohen-Adad, S. Chapuisat, J. Doyon, S. Rossignol, J. M. Lina, H. Benali, and F. Lesage, “Activation detection in diffuse optical imaging by means of the general linear model,” Med. Image Anal. 11(6), 616–629 (2007). [CrossRef]
8. M. S. Hassanpour, B. R. White, A. T. Eggebrecht, S. L. Ferradal, A. Z. Snyder, and J. P. Culver, “Statistical analysis of high density diffuse optical tomography,” NeuroImage 85(1), 104–116 (2014). [CrossRef]
9. A. Aarabi, V. Osharina, and F. Wallois, “Effect of confounding variables on hemodynamic response function estimation using averaging and deconvolution analysis: An event-related NIRS study,” NeuroImage 155, 25–49 (2017). [CrossRef]
10. A. K. Seghouane and A. Shah, “HRF estimation in fMRI data with an unknown drift matrix by iterative minimization of the Kullback-Leibler divergence,” IEEE Trans. Med. Imaging 31(2), 192–206 (2012). [CrossRef]
11. W. Chen, “Hemodynamic Response Function Modeling,” (University of North Carolina, 2012).
12. T. Zhang, F. Li, L. Beckes, and J. A. Coan, “A semi-parametric model of the hemodynamic response for multi-subject fMRI data,” NeuroImage 75, 136–145 (2013). [CrossRef]
13. M. A. Lindquist and T. D. Wager, “Validity and power in hemodynamic response modeling: a comparison study and a new approach,” Hum. Brain Mapp. 28(8), 764–784 (2007). [CrossRef]
14. F. Pedregosa, M. Eickenberg, P. Ciuciu, B. Thirion, and A. Gramfort, “Data-driven HRF estimation for encoding and decoding models,” NeuroImage 104, 209–220 (2015). [CrossRef]
15. G. Marrelec, P. Ciuciu, M. Pelegrini-Issac, and H. Benali, “Estimation of the hemodynamic response in event-related functional MRI: Bayesian networks as a framework for efficient Bayesian modeling and inference,” IEEE Trans. Med. Imaging 23(8), 959–967 (2004). [CrossRef]
16. K. J. Friston, O. Josephs, G. Rees, and R. Turner, “Nonlinear event-related responses in fMRI,” Magn. Reson. Med. 39(1), 41–52 (1998). [CrossRef]
17. D. A. M. J. H. B. Mapping, “Optimal experimental design for event-related fMRI,” Hum Brain Mapping8 (1999).
18. R. Casanova, S. Ryali, J. Serences, L. Yang, R. Kraft, P. J. Laurienti, and J. A. Maldjian, “The impact of temporal regularization on estimates of the BOLD hemodynamic response function: a comparative analysis,” NeuroImage 40(4), 1606–1618 (2008). [CrossRef]
19. J. J. Shynk, Probability, Random Variables and Random Processes Theory and Signal Processing Applications (John Wiley & Sons, Inc., 2013), p. 291.
20. L. H. Nguyen and K.-S. Hong, “Investigation of the Hemodynamic Response in Near Infrared Spectroscopy Data Analysis,” 2010 Second International Conference on Knowledge and Systems Engineering, 28–32 (2010). [CrossRef]
21. N. M. Gregg, B. R. White, B. W. Zeff, A. J. Berger, and J. P. Culver, “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front. Neuroenerg. 2, 14 (2010).
22. N. Bazargani and A. Nosratinia, “Joint maximum likelihood estimation of activation and hemodynamic Response Function for fMRI,” Med. Image Anal. 18(5), 711–724 (2014). [CrossRef]
23. S. Kohno, I. Miyai, A. Seiyama, I. Oda, A. Ishikawa, S. Tsuneishi, T. Amita, and K. Shimizu, “Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis,” J. Biomed. Opt. 12(6), 062111 (2007). [CrossRef]
24. J. Markham, B. R. White, B. W. Zeff, and J. P. Culver, “Blind identification of evoked human brain activity with independent component analysis of optical data,” Hum. Brain Mapp. 30(8), 2382–2392 (2009). [CrossRef]
25. C. Eckart and G. Young, “The approximation of one matrix by another of lower rank,” Psychometrika 1(3), 211–218 (1936). [CrossRef]
26. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. A. Sato, “Extended hierarchical Bayesian diffuse optical tomography for removing scalp artifact,” Biomed. Opt. Express 4(11), 2411–2432 (2013). [CrossRef]
27. B. Wang, T. Pan, Y. Zhang, D. Liu, J. Jiang, H. Zhao, and F. Gao, “A Kalman-based tomographic scheme for directly reconstructing activation levels of brain function,” Opt. Express 27(3), 3229–3246 (2019). [CrossRef]
28. Q. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plucker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef]
29. V. Fonov, A. C. Evans, K. Botteron, C. R. Almli, R. C. McKinstry, and D. L. Collins, “Unbiased average age-appropriate atlases for pediatric studies,” NeuroImage 54(1), 313–327 (2011). [CrossRef]
30. A. T. Eggebrecht, B. R. White, S. L. Ferradal, C. Chen, Y. Zhan, A. Z. Snyder, H. Dehghani, and J. P. Culver, “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” NeuroImage 61(4), 1120–1128 (2012). [CrossRef]
31. V. Bonomini, L. Zucchelli, R. Re, F. Ieva, L. Spinelli, D. Contini, A. Paganoni, and A. Torricelli, “Linear regression models and k-means clustering for statistical analysis of fNIRS data,” Biomed. Opt. Express 6(2), 615–630 (2015). [CrossRef]
32. G. H. Glover, “Deconvolution of impulse response in event-related BOLD fMRI,” NeuroImage 9(4), 416–429 (1999). [CrossRef]
33. Y. Zhang, D. H. Brooks, and D. A. Boas, “A haemodynamic response function model in spatio-temporal diffuse optical tomography,” Phys. Med. Biol. 50(19), 4625–4644 (2005). [CrossRef]
34. F. Scarpa, S. Brigadoi, S. Cutini, P. Scatturin, M. Zorzi, R. Dell’acqua, and G. Sparacino, “A reference-channel based methodology to improve estimation of event-related hemodynamic response from fNIRS measurements,” NeuroImage 72, 106–119 (2013). [CrossRef]
35. T. Sato, I. Nambu, K. Takeda, T. Aihara, O. Yamashita, Y. Isogaya, Y. Inoue, Y. Otaka, Y. Wada, M. Kawato, M. A. Sato, and R. Osu, “Reduction of global interference of scalp-hemodynamics in functional near-infrared spectroscopy using short distance probes,” NeuroImage 141, 120–132 (2016). [CrossRef]
36. C. Habermehl, C. Schmitz, S. P. Koch, J. Mehnert, and J. Steinbrink, “Investigating hemodynamics in scalp and brain using high-resolution diffuse optical tomography in humans,” BSu2A.2 (2012).
37. J. R. Goodwin, C. R. Gaudet, and A. J. Berger, “Short-channel functional near-infrared spectroscopy regressions improve when source-detector separation is reduced,” Neurophotonics 1(1), 015002 (2014). [CrossRef]
38. 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–1882 (2005). [CrossRef]
39. B. Vojnovic, “Advanced time-correlated single photon counting techniques,” J. Microsc. 222(1), 65–66 (2006). [CrossRef]
40. D. Qin, Z. Ma, F. Gao, and H. Zhao, “Determination of optical properties in turbid medium based on time-resolved determination,” Fifth International Conference on Photonics and Imaging in Biology and Medicine (2007).
41. D. Liu, B. Wang, T. Pan, J. Li, Z. Qin, L. Zhang, Z. Zhou, and F. Gao, “Toward quantitative near infrared brain functional imaging: lock-in photon counting instrumentation combined with tomographic reconstruction,” IEEE Access 7, 86829–86842 (2019). [CrossRef]
42. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection: methodology and phantom validation,” Phys. Med. Biol. 49(6), 1055–1078 (2004). [CrossRef]
43. C. Chen, F. Tian, H. Liu, and J. Huang, “Diffuse optical tomography enhanced by clustered sparsity for functional brain imaging,” IEEE Trans. Med. Imaging 33(12), 2323–2331 (2014). [CrossRef]
44. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse recovery methods hold promise for diffuse optical tomographic image reconstruction,” Ieee J Sel Top Quant 20(2), 74–82(2014). [CrossRef]
45. C. Habermehl, J. Steinbrink, K. R. Muller, and S. Haufe, “Optimizing the regularization for image reconstruction of cerebral diffuse optical tomography,” J. Biomed. Opt. 19(9), 096006 (2014). [CrossRef]
46. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]
47. G. H. Golub and C. F. V. Loan, Matrix Computations (The Johns Hopkins University Press, 2013).
48. R. A. Horn and C. R. Johnson, Matrix Analysis, second ed. (Cambridge University Press, 2012).
