Joint direct estimation of hemodynamic response function and activation level in brain functional high density diffuse optical tomography

: High density diﬀuse optical tomography has become increasingly important to detect underlying neuronal activities. Conventional methods ﬁrst estimate the time courses of the changes in the absorption coeﬃcients 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 ﬁnal 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 ﬂuxes. It is assumed that all the voxels in the same activated brain region share the same HRF but diﬀer in the activation levels, and no prior information is imposed on the speciﬁc 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.


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][4][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 (µ a ) at every time point are separately reconstructed and a time course for every voxel is formed. Then the changes in µ 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 µ 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 µ a changes instead of being directly estimated from the boundary fluxes. The error in the separate computation of the µ 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 µ 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.

Methods
The inverse problems in DOT at all the time points (k=1, 2, . . . . . . , K) can be stacked together as follows: where M = [m 1 m 2 · · · m K ] ∈ R M×K with m k being the boundary flux changes at the k-th time point and M being the total number of the measurements; X = [x 1 x 2 · · · x K ] ∈ R N×K with x k being the µ a changes at the k-th time point and N being the total number of the voxels; J ∈ R M×N is the Jacobian matrix to describe the sensitivity of m k to x k . According to the GLM, X can be expressed in the following form: where β = [β 1 β 2 · · · β N ] T with β n denoting the activation level at the n-th voxel; (·) T denotes the transpose of the vector or matrix; S is the stimulus convolution matrix; Its entry S(i, j) will be 1 if a stimulus is on at time k = i − j, k ≥ 0, and 0 otherwise. It is assumed that the neuronal cells in the joint brain regions have the similar functions and share the same HRF h of µ a in response to the same stimulus. That is to say, h = h n for n=1, 2, . . . . . . , N. α = [α 1 α 2 · · · α N ] T with α n denoting the amplitude of the nuisance interference signal f at the n-th voxel; E is the noise term.
Equation (1) can be further reformulated as follows by substituting Eq. (2) into it where E 1 = JE is the noise term. 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 (f) and hemodynamic signals in the brain (h) are usually assumed to be independent from each other [20,21]. f denotes the periodic changes of µ 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. h and f can be further inferred to be orthogonal each other.
Suppose that r(·) denote the rank of a vector or matrix. r(β), r(h), r(α) and r(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(Jβ(Sh) T ) = 1 and r(Jα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(M) ≤ 2. On the other hand, because M has at least two independent components h and f, r(M) ≥ 2 [23,24]. Therefore, we obtain r(M) = 2.
In the singular value decomposition of 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(M) = 2, the problem in Eq. (3) is essentially to find a rank-2 approximation to the measurement M. According to the Eckart and Young theorem, a rank-2 approximation can be obtained by the truncated singular value decomposition method [25]. Specifically, M can be approximated by where ξ i , u i and 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. M is more sensitive to the scalp interference signals f than the cerebral signals h induced by the brain activation [26]. Thus, f has bigger weight than h. Therefore, the first term and second term in Eq. (4) respectively account for Jαf T and Jβ(Sh) T . Accordingly, the second term can be expressed as follows: By further decomposing Eq. (5), it is easy to obtain Jβ = u 2 and Sh = ξ 2 v 2 . By solving separately them, β and 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 h is a smooth curve. This a prior information is helpful to regularize h. The equation Sh = ξ 2 v 2 is converted into the following optimization problem by the commonly used Tikhonov regularization method [22]: where λ>0 is a positive regularization parameter tuning the balance between the fidelity and the regularization terms. The second order difference matrix C defined below is used to impose the smoothness constraint

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 h est and h tr , respectively. Then, the root mean square error (RMSE) is defined below [22] Besides, to quantitatively compare the qualities of the activation-level maps estimated by different methods, several metrics are adopted [27]. First, the RMSE in this case is defined as follows where β est and β tr are the estimated and true activation-level maps (both normalized by their own maximum values). Second, the contrast-to-noise ratio (CNR) is selected to evaluate how the activated region can be recognized from the background. The voxels with β bigger than half of the maximum β value are selected to compose the activated brain region. The other voxels compose the background. Then CNR is defined as follows where β atd and β bk are the vectors consisting of the β for the voxels in the activated region and background, respectively; w = |β atd |/(|β atd | + |β bk |) with | · | representing the number of the elements in the vector. Area ratio (AR) is adopted to assess the ratio of the area of the estimated activated region to the true one. At last, the duration is used to compare the computation speed. In general, smaller RMSE, bigger CNR and AR closer to 1 represent higher quality images. Shorter duration means faster computation speed.

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 mm, 77 mm < Y < 85 mm and 80 mm < Z < 100 mm consist the activated brain region, as shown in Fig. 1 and Fig. 2.  The HRF of HbO concentration is modeled as the linear combination of two gamma functions [31,32]: The green color denotes the activated brain region. where k denotes the discrete time points; Γ n (k, τ j , ρ j ) = 1 ; ψ 1 regulates the amplitude; ψ 2 determinates the undershoot; τ 1 and τ 2 regulate the h HbO (k) shape; ρ 1 and ρ 2 tune the scale. These specific parameters are set as follows: The HRF for the HbR h HbR (k) is inverted and with a maximum set at -1/3 with respect to h HbO (k). Besides, h HbR (k) delays 2 s relative to h HbO (k). In summary, h HbR (k) = − 1 3 h HbO (k − 2f s ) with f s denoting the sampling rate [33]. Their waveforms are shown in Fig. 3(a).
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 µ a changes at the wavelength of 760 nm and 830 nm are calculated and shown in Fig. 3(c). Besides to the changes in µ 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 µ 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][35][36][37]. For the local case, the interferences are assumed to happen only in the scalp in the domain 24mm<X<50 mm, 50 mm < Y < 100 mm and 40 mm < Z < 75 mm.
The NN-1 channels with a short source-detector separation are only sensitive to the µ 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 µ 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 Butterworthtype 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.

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 µ a in the activated brain region. These observations confirm the conclusions in Section 2.
The time courses of changes in µ 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.   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 RMSE β , 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.

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 µ a in gray matter. Then the changes of µ 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.  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.

Phantom experiments
The performances of the proposed method are further investigated based on dynamic phantom experiments [27].

Experimental setup
As shown in Fig. 10, a cuboid-shaped daicel phantom of size 100 mm × 120 mm × 30 mm is used to simulate a human head. At the wavelength of 785 nm, its optical properties are about µ a = 0.0041/mm, µ 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 (µ a = 0.0260/mm, µ s = 1.0134/mm) and background solution (µ a = 0.0173/mm, µ 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, µ a in the cylinder hole will first increase and then decrease back to the initial value. 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.

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 µ a changes (δµ 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. RMSE h of JDE is smaller than that of SFIR, as shown in the first row in Table 8. 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 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 RMSE β , bigger CNR and AR closer to 1. The last row of Table 8 shows that the JDE has a shorter duration than SFIR. The metrics of the normalized activation-level maps are compared in Table 8 where the JDE obtains smaller RMSE β , bigger CNR and AR closer to 1. The last row of Table 8 shows that the JDE has a shorter duration than SFIR.  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 RMSE β , bigger CNR and AR closer to 1. The last row of Table 8 shows that the JDE has a shorter duration than SFIR.

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 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 u 1 and 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 µ 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. RMSE β , CNR and AR are respectively 0.52, 17.90 and 0.91.

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 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 1 u and 1 v . 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 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

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 µ 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