Translator Disclaimer
1 November 2009 Monte Carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy
Author Affiliations +
Abstract
The performance of near-infrared spectroscopy is sometimes degraded by the systemic physiological interference in the extracerebral layer. There is some systemic interference, which is highly correlated with the functional response evoked by a task execution. This kind of interference is difficult to remove by using ordinary techniques. A multidistance measurement method is one of the possible solutions for this problem. The multidistance measurement method requires estimation parameters derived from partial pathlength values of tissue layers to calculate an absorption coefficient change from a temporal absorbance change. Because partial path lengths are difficult to obtain, experimentally, we estimated them by a Monte Carlo simulation based on a five-layered slab model of a human adult head. Model parameters such as thickness and the transport scattering coefficient of each layer depend on a subject and a measurement position; thus, we assumed that these parameters obey normal distributions around standard parameter values. We determined the estimation parameters that provide a good separation performance in average for the model parameter distribution. The obtained weighting is robust to model parameter deviation and provides smaller errors on average compared to the parameters, which are determined without considering parameter distribution.

1.

Introduction

Near-infrared spectroscopy (NIRS) has several advantages compared to other functional measurement methods, such as functional magnetic resonance imaging (fMRI), positron emission tomography, and electroencephalography (EEG). The advantages include good temporal resolution and measurement of both oxygenated (HbO2) and deoxygenated (HbR) hemoglobin, portability, low restraints of subjects, and low cost of the measurement equipment. The disadvantages are low spatial resolution and potentially unstable measurement because of the instability of an optical probe contact that might be induced by body motion and hair absorption. The systemic physiological interference is an another difficult problem in NIRS measurement.

The well-known systemic physiological interference is caused by heart beat and respiration. The suppression of this kind of systemic interference is relatively easy because the correlation between the interference and the functional response of brain functional activity is usually low. For example, systemic interference caused by heart beat can be effectively removed using a low pass filtering technique. However, low pass filtering may not be appropriate for respiration because the frequency spectra of respiration and functional activity may overlap each other. Even in this case, it can be suppressed by using the low correlation property between the interference (respiration) and brain functional activity. Zhang solved this problem by using an adaptive filtering technique.1, 2 Saager and Berger applied the least-squares fit technique to this problem, and gave a detailed analysis of its performance.3 Both methods are based on the multidistance measurement and have an advantage that they can remove not only the interference in the superficial layer (scalp and skull layers), but also the interference included in the cerebral layer, as far as the interference is low correlated with the brain functional signal.

There is an another form of interference, which is observed synchronously with the heart rate increase during execution of a task (for example, finger opposition task).4, 5 This kind of interference is highly correlated with the functional response itself. Thus, it is difficult to remove it by the methods using the low correlation property between them. We consider the removal of this kind of interference in this paper.

Two approaches have been proposed to isolate and eliminate this kind of interference. The first approach is to use the spatial locality difference between the brain functional activity and the systemic interference. It is well known through fMRI studies that the brain functional activity is highly localized. Conversely, systemic interference is generally global. Thus, this difference can be used to differentiate them. Pioneering work of this approach was reported by Franceschini 4 They assumed that the NIRS signal measured at an area where no functional activity is expected should be similar to the systemic interference, and isolated the functional response of the activated region by subtracting this signal from the measured signal. This approach was refined as a method that uses principal component analysis.5, 6 A similar approach based on independent component analysis instead of PCA has also been proposed.7 This approach has an advantage in that a conventional continuous-wave NIRS system is available, and the method is easy to implement. However, it has a tendency to decrease the amplitude of the functional response in the activated regions and to propagate noise from noisy channels to all other channels.5 Another disadvantage is that it requires a large-scale NIRS system equipping many measurement channels because it is based on the spatial uniformity of the systemic interference.

The second approach relies on the multidistance measurement. If we assume a multi layered slab model as a human adult head, the temporal absorbance change can be approximated as a linear sum of the products of a partial pathlength and an absorption coefficient change of each layer. 8, 9 Because the partial pathlength depends on the distance between a light source and a detector, absorption coefficient change at each layer can be computed easily if several detectors with different source-detector separations are used and partial pathlength parameters of each source-detector pair are known. 8, 10 Thus, if such an interference occurs mainly in the superficial layer, multidistance measurement method can separate it from the functional response.

The problem with the second approach is that it is very difficult to determine, experimentally, the partial pathlength of each layer. Fabbri 10 pointed out that the ratio of the partial pathlength of the superficial layer at different detectors can be calculated as the ratio of the absorbance changes, if the hemoglobin changes in the cerebral layer are negligible. This approach is effective. However, requirement of static hemodynamics in the cerebral layer may be a strong constraint in some situations.

A partial pathlength of each layer can be predicted computationally by Monte Carlo simulation. Monte Carlo simulation requires a lot of physiological parameters of a head model, such as the thickness and transport scattering coefficient of each tissue. Tissue thickness can be measured by MRI, however, doing so may spoil the advantages of low restraint of subjects and low cost provided by NIRS. For example, although NIRS offers advantages for brain functional measurement of infants and children because of its low restraint property, its applicability would be limited if MRI measurement is required. The transport scattering coefficient of each layer is much more difficult to measure.

Leung 11 have provided a theoretical analysis of the second approach by the finite element method of the diffusion equation based on assumed standard physiological parameters. They showed how exactly the hemoglobin change at gray matter can be estimated using this approach and gave the optimum detector arrangement. In this paper, deviations of blood volume, oxygen saturation, and thickness of a superficial layer are considered, and a range of data were simulated corresponding to different thicknesses and physiological conditions of each layer. The optimum regression parameters for brain hemodynamic change estimation were calculated directly from the simulated data without using the modified Beer-Lambert law. They used the partial least-squares method to obtain these parameters. In their paper, however, only the thickness parameter of a superficial layer was considered as a deviation, and the transport scattering parameter of all layers and thickness parameter of the other layers were not considered.

Here, we analyzed the multidistance measurement method to separate the brain functional response from the systemic interference in a superficial layer caused by task execution. Other interferences, evoked by heart beat, respiration, etc., were assumed to be removed beforehand by preprocessing. We used the modified Beer-Lambert law, and partial optical path lengths were estimated by Monte Carlo simulation. We supposed a population of subjects and assumed that physiological parameters, the transport scattering parameters, and thicknesses of layers of each subject were distributed around the assumed standard values. An estimation method giving a good average separation performance for the population was derived by multilinear regression.

2.

Theory

Our study used two light detectors ( d1 and d2 ) and two wavelengths ( λ1 and λ2 ) of light. A five-layer slab model consisting of scalp, skull, cerebrospinal fluid (CSF), gray matter, and white matter was used as a human adult head model.

2.1.

NIRS Observation Model Based on Partial Pathlength

The temporal absorbance change at the detector di (i=1,2) of wavelength λj (j=1,2) is represented as follows:8, 9

Eq. 1

ΔAdi,λj=ldi,λjscΔμa,λjsc+ldi,λjskΔμa,λjsk+ldi,λjcsfΔμa,λjcsf+ldi,λjgmΔμa,λjgm+ldi,λjwmΔμa,λjwm,
where ldi,λjsc , ldi,λjsk , ldi,λjcsf , ldi,λjgm , and ldi,λjwm are partial path lengths of scalp, skull, CSF, gray matter, and white matter, respectively, and Δμa,λjsc , Δμa,λjsk , Δμa,λjcsf , Δμa,λjgm , and Δμa,λjwm are absorption coefficient changes of each layer. We assume that absorption coefficient change in CSF is relatively small (Δμa,λjcsf0) because the concentrations of HbO2 and HbR are small there. We also assume ldi,λjwm0 because the result of Monte Carlo simulation shows that the partial path length of white matter is very small. Finally, we assume that absorption coefficient changes in scalp and skull are similar (Δμa,λjscΔμa,λjsk) . Equation 1 thus becomes

Eq. 2

ΔAdi,λj=ldi,λjspΔμa,λjsp+ldi,λjgmΔμa,λjgm,
where ldi,λjsp and Δμa,λjsp are partial pathlength and absorption coefficient change, respectively, in superficial layer (scalp and skull)

Eq. 3

ldi,λjsp=ldi,λjsc+ldi,λjsk,

Eq. 4

Δμa,λjsp=Δμa,λjsc=Δμa,λjsk.
If we give a vector representation of the temporal absorbance change and absorption coefficient change, then we have

Eq. 5

a=(ΔAd1,λ1,ΔAd2,λ1,ΔAd1,λ2,ΔAd2,λ2)T,

Eq. 6

μ=(Δμa,λ1sp,Δμa,λ1gm,Δμa,λ2sp,Δμa,λ2gm)T.
Equation 2 is summarized as follows:

Eq. 7

a=Lμ,
where

Eq. 8

L=(ld1,λ1spld1,λ1gm00ld2,λ1spld2,λ1gm0000ld1,λ2spld1,λ2gm00ld2,λ2spld2,λ2gm).
We call L a partial pathlength matrix.

Next we give a vector representation of HbO2 and HbR concentration change in superficial and gray matter layers

Eq. 9

x=(ΔHbO2sp,ΔHbO2gm,ΔHbRsp,ΔHbRgm)T.
Hemoglobin concentration change x and absorption coefficient change μ have the following relationship:

Eq. 10

μ=Ex.
E is a molar absorption coefficient matrix of HbO2 and HbR

Eq. 11

E=(εHbO2,λ10εHbR,λ100εHbO2,λ10εHbR,λ1εHbO2,λ20εHbR,λ200εHbO2,λ20εHbR,λ2).
From Eqs. 7, 10, we derive the following.

Eq. 12

a=LEx.

2.2.

Optimizing the Estimation of ΔHbO2 and ΔHbR

Hemoglobin concentration change can be estimated from the measured temporal absorbance change a and Eq. 12 if we know E and L . E can be obtained from prior studies.12 However, partial pathlength matrix L is determined by physiological parameter values, which depend on a subject and a measurement position, and it is difficult to obtain experimentally. We describe in this paper a separation method that provides a good average separation performance for a variety of subjects and measurement positions.

The various partial pathlength matrices are written as L={L1,L2,,LK} corresponding to a variety of subjects and measurement positions. L indicates that various values of the temporal absorbance change ak=LkEx are observed even if the hemoglobin concentration change x is fixed. We introduce an inverse transformation matrix E1W , which estimates hemoglobin concentration change from temporal absorbance change values. W estimates absorption coefficient change from the temporal absorbance change, and E1 transforms absorption coefficient change to hemoglobin concentration change. Here we assume that the absorption coefficient changes of two wavelengths are estimated independently. Thus, W=diag(W1,W2) , where W1 and W2 are estimation matrices for λ1 and λ2 , respectively. Because W includes coefficients to estimate absorption coefficient changes, we call W a measurement weighting.

When we use a measurement weighting W , the hemoglobin concentration change when L=Lk is estimated as follows:

Eq. 13

x̂k(t)=E1Wak(t),

Eq. 14

=E1WLkEx(t).
The estimation error is

Eq. 15

x̂k(t)x(t)=(E1WLkEI)x(t).
If we assume a model time course of hemodynamic change, then mean-squared error of estimated hemoglobin changes is given as follows:

Eq. 16

Jk=1T0T[x̂k(t)x(t)]T[x̂k(t)x(t)]dt,

Eq. 17

=1T0Tx(t)T(E1WLkEI)T(E1WLkEI)x(t)dt,

Eq. 18

=tr{(E1WLkEI)T(E1WLkEI)Σx},
where

Eq. 19

Σx=1T0Tx(t)xT(t)dt.
Jk is a mean-squared error of estimated hemoglobin changes when the partial path-length matrix is Lk . Thus, the mean-squared error on average for all partial path-length matrices in L is given as follows:

Eq. 20

J=1Kk=1KJk.

The optimum weighting W=diag(W1,W2) to minimize J and its derivation detail is given in the Appendix. The optimum weighting [Eq. 43] gives a good average separation performance for a partial path-length matrix distribution L on average. Thus, if we apply this weighting to the multidistance measurement data, it is expected to give a better separation performance than the weighting determined without considering a parameter distribution. This will be illustrated later in Fig. 5.

Fig. 5

Mean-squared errors of ΔHbOgm and ΔHbRgm estimated by Wopt (filled bars) and Lstd1 (unfilled bars). Standard deviations of normal distributions of model parameters were set to be 10, 20, and 30% of standard parameter values shown in Table 1. (a) and (c) MSEs of ΔHbOgm , (b) and (d) MSEs of ΔHbRgm , and (a) and (b) standard artifact case. Changes in the superficial layer is 110th of changes in the gray-matter layer. (c) and (d) Large artifact case. Changes in the superficial layer is 15 of changes in the gray-matter layer.

064025_1_037906jbo5.jpg

2.3.

Random Generation of Partial Path-length Matrices

Because a partial path-length matrix L is a block diagonal matrix, we write L=diag(Λ1,Λ2) . Λj corresponds to a measurement of λj light. The value of a partial path-length depends on the thickness and transport scattering coefficient of each layer of a head. Thus, we use Λj(uj) to emphasize this. The set of parameters uj are as follows:

Eq. 21

uj=(u1j,u2j,u3j,u4j,u5j,u6j,u7j,u8j,u9j)T

Eq. 22

=(μs,λjsc,μs,λjsk,μs,λjcsf,μs,λjgm,μs,λjwm,tsc,tsk,tcsf,tgm)T.
The thickness of white matter is not included in uj because our Monte Carlo simulation shows that few photons can reach white matter. The value of a partial pathlength also depends on the baseline concentrations of HbO2 and HbR; however, we assume here that the baseline concentration is fixed at its standard value to simplify the analysis.

A set of partial path-length matrices can be obtained if we execute a Monte Carlo simulation based on various multilayered model parameters. This simulation, however, requires significant computation time; thus, execution of many sets of parameters is difficult. Instead, we assumed that Λj(u0j+Δuj) can be approximated in a linear form if Δuj is small

Eq. 23

Λj(u0j+Δuj)=Λj(u0j)+m=19|Λjum|u0jΔumj.

Many partial path-length matrices can be generated based on Eq. 23. The first-order coefficients of Eq. 23 indicate how much the partial path length changes when thicknesses and transport scattering coefficients of layers change; we term this the sensitivity coefficient of the partial pathlength. The approximate value of the sensitivity coefficient can be obtained if we investigate the partial pathlength change produced by a small increase in the parameter value um through simulations.

We computed the partial-pathlength and sensitivity coefficient of standard parameters um,0j by Monte Carlo simulation. We assumed that a parameter umj is a random variable and obeys a normal distribution N(um,0j,σmj) . If we determine a parameter value umj by random sampling of N(um,0j,σmj) , then many partial pathlength matrices can be generated by Eq. 23.

3.

Method

3.1.

Monte Carlo Simulation

The partial path-length and sensitivity coefficient based on standard parameter values were computed by a Monte Carlo simulation. The simulation program, tMCimg, provided in public by the Photon Migration Imaging Laboratory at Massachusetts General Hospital was used. A five-layer model of an adult head consisting of scalp, skull, CSF, gray matter, and white matter was used (Fig. 1 ). The size of simulated tissues was 100×100×50mm . The wavelengths of the light source were 800 and 840nm . The source-detector separation was assumed from 10to40mm every 5mm .

Fig. 1

Geometry of a five-layer adult head model for Monte Carlo simulation.

064025_1_037906jbo1.jpg

The standard parameters used in the simulation are given in Table 1 . Thickness and optical properties ( μa and μs ) for near-infrared light of 800nm for each layer were chosen from reported data.13 HbO2 and HbR concentrations were calculated from the absorption coefficients at 800nm with the assumption that oxygen saturation is 70%. The absorption coefficient calculated from water absorption coefficients given in the paper14 was also used in this case.15 These values were used to calculate the absorption coefficients at 840nm . Because it is known that the transport scattering coefficient decreases if the wavelength of infrared light increases,11, 16 transport scattering coefficients at 840nm were set to be 95% of those at 800nm .

Table 1

Hemodynamic parameters, thickness, and optical properties for each layer of an adult head model.

Tissue typeBaselineconcentrationThickness(mm)Absorptioncoefficient μa 800nm∕840nm (mm−1) Transport scatteringcoefficient μs′ 800nm∕840nm (mm−1)
HbO2 (mM)HbR(mM)
Scalp0.0640.0273 0.0180.021 1.901.81
Skull0.0570.0247 0.0160.019 1.601.52
CSF0.0140.0062 0.0040.005 0.240.23
Gray matter0.1280.0554 0.0360.042 2.202.09
White matter0.0500.02134 0.0140.017 9.108.65

3.2.

Model Time Course of Hemodynamic Change

It is very difficult to know real values x=(ΔHbOsp,ΔHbO2gm,ΔHbRsp,ΔHbRgm)T in superficial and gray-matter layers. Thus, we use a simple model of hemodynamic change time course. In this model, ΔHbO2gm and ΔHbRgm are functional responses of the executed task, and ΔHbOsp and ΔHbRsp are systemic interference correlated with the task execution.

The hemodynamic responses in both layers were defined as the convolution of the stimulation s(t) , [ s(t)=0 for rest period, and 1 for stimulation] and a prototypical hemodynamic impulse response h(t) ,1, 17 in our model

Eq. 24

s(t)={1iftrest0iftstimulation},

Eq. 25

h(t)=tbexp(bDt).
The parameters b and D in the definition of h(t) determine the width and the delay time to the peak point of the impulse response, respectively. The delay time to the peak point for both ΔHbOgm and ΔHbRgm in the gray-matter layer were set D=5s , and these parameters for the superficial layer were set D=3s for both ΔHbOsp and ΔHbRsp because the recent findings18 show that ΔHbOgm and ΔHbRgm have a similar time course with a very little delay and that the functional response is slower than the systemic response. The parameter concerning the impulse response width was set at b=8 for all responses. Figure 2a shows the stimulation s(t) (0t40) , and Fig. 2b shows the hemodynamic impulse response of each layer. The peak values of all responses are normalized to be unity.

Fig. 2

Model time course of hemodynamic change: (a) Stimulation paradigm, (b) normalized hemodynamic impulse response, and (c) evoked hemodynamic change.

064025_1_037906jbo2.jpg

The evoked hemodynamic response x(t) was the convolution of s(t) and h(t)

Eq. 26

x(t)=C[h(t)s(t)],
where C is a scaling parameter. Figure 2c shows the simulated evoked hemodynamic response of each layer. Peak values of ΔHbOgm and ΔHbRgm were set to be 0.017 and 0.005mM , respectively. These values were about 13% and 9% of the assumed baseline conditions.19 On the other hand, peak values of the superficial layer were set to be 110 of the gray-matter layer based on the following findings: The temporal absorbance changes of the systemic interference of a finger-tapping task and functional response of brain functional activity are considered to be comparable4 and Monte Carlo simulation shows that the partial path length of a superficial layer is more than 10 times longer than that of a gray-matter layer. Thus, we assumed that hemoglobin concentration changes of a superficial layer are 110 of those of gray-matter layer.

The covariance matrix Σx [Eq. 19] was calculated based on these simulated hemodynamic responses and used for the optimization process.

4.

Results and Discussion

4.1.

Partial Path Length and Its Sensitivity Coefficient

We computed partial path lengths of superficial and gray-matter layers at various source-detector separations by Monte Carlo simulation (Fig. 3 ). Solid lines and dashed lines in Fig. 3 correspond to the light sources of wavelengths λ1 and λ2 , respectively. Figure 3 shows that partial path lengths of different wavelength have similar values, especially for gray matter. It also shows that the partial path length of the superficial layer is more than 10 times greater than that of gray matter. This difference is pronounced when the source-detector separation is small, which means that blood flow change in a superficial layer may have a strong influence on the measurement.

Fig. 3

Partial path length of superficial and gray-matter layer based on the standard model parameters. Solid lines indicate the partial pathlength of λ1 , and dashed lines indicate the partial pathlength of λ2 .

064025_1_037906jbo3.jpg

Figure 4 shows relative (percent) change of partial pathlength of gray-matter and superficial layers, that is, the ratio of partial path-length change to its baseline value, where model parameters were increased by 10% from the standard values. We showed these values instead of sensitivity coefficients because they may be more comprehensible. Figure 4 indicates that the thickness parameter of scalp and skull layers has a significant influence on the partial path length of the gray-matter and superficial layers. The transport scattering coefficient of the skull layer has also a comparable influence on the partial pathlength of the gray-matter layer. On the other hand, partial path lengths of both layers are insensitive to all other parameters.

Fig. 4

Relative (percent) change of partial path length in superficial and gray-matter layers, where model parameters were increased by 10% from the standard values: (a) Cases in which the transport scattering coefficient of each layer was increased and (b) cases in which the thickness of each layer was increased. X -axis corresponds to the source-detector separation.

064025_1_037906jbo4.jpg

The validity of the first-order approximation of the partial path-length matrix [Eq. 23] was checked by another Monte Carlo simulation. The simulation result (data not shown) shows Eq. 23 holds approximately in the range of ±30% of the transport scattering coefficient change and ±1mm of the layer thickness change.

4.2.

Optimum Weighting and Its Performance

To obtain a measurement weighting that provides a good average performance for the subject population, we generated a set L of partial path-length matrices based on its linear approximation model [Eq. 23] and the assumption that model parameters are distributed normally around the standard parameter values. The source-detector separations of two detectors were set to be 20 and 30mm . The optimum weighting Wopt for L was derived by Eq. 43. Standard deviations of normal distributions of model parameters were set to be 10, 20, and 30% of standard parameter values shown in Table 1. Two thousand partial path-length matrices were randomly generated for each case, and half of them were used as training data to compute Wopt and the other half as test data to evaluate the performance of the obtained Wopt .

When we use a simulated evoked hemodynamic response of Fig. 2, the mean-squared error (MSE) of ΔHbOgm and ΔHbRgm estimated by two measurement weightings, Wopt and Lstd1 are shown in Figs. 5a and 5b. Wopt is the optimum weighting obtained by the training data, and its MSEs are shown by filled bars. Lstd1 is the inverse of the standard partial path-length matrix Lstd , where Lstd is a partial path-length matrix determined by the standard model parameters. This means that Lstd1 is a weighting not to consider the model parameter distribution. Its MSEs are shown by unfilled bars. Figures 5a and 5b show that MSE of both ΔHbOgm and ΔHbRgm increases almost linearly when deviations of model parameters increases. It also shows that MSEs of the optimum weighting are much less than those obtained by Lstd1 , especially when the deviation of model parameter distribution is large. Thus, the proposed optimization method is effective to minimize the estimation errors on average.

Figures 5c and 5d show MSEs of the case where the hemoglobin changes in the superficial layer is larger than the case of Fig. 2. In this case, the peak values of ΔHbOsp and ΔHbRsp were set 15th of the gray matter layer values. Figures 5c and 5d show that MSEs of the optimum weighting slightly increase in this case; however, MSEs of Lstd1 are significantly larger than the small artifact case (Fig. 2 case). This result shows that the optimum weighting can give a better and robust estimation of brain hemodynamic response on average.

Figure 6 shows how much the estimated hemoglobin changes deviate from the assumed original changes (Fig. 2) when model parameter changes. The standard deviation of model parameter distribution was assumed to be 20% of the standard values. Ten partial path-length matrices were randomly selected from the generated test data, and two measurement weightings ( Wopt and Lstd1 ) were used to estimate the hemoglobin changes based on Eq. 14. The estimated hemoglobin changes are light red (HbO) and blue (HbR) solid lines. The assumed original hemoglobin changes are also shown in the figures by thick solid lines. The estimation result by Wopt and Lstd1 are shown in Figs. 6a and 6b, respectively, which show that the estimated changes by Lstd1 vary widely compared to those by Wopt .

Fig. 6

(a, b) Hemoglobin changes estimated by two weightings ( Wopt and Lstd1 ) and dual-detector NIRS algorithm. (c) Hemoglobin changes estimated by an ordinary single detector NIRS algorithm. (d–f) Normalized hemoglobin changes of (a–c). Ten different partial path-length matrices were used in the simulation, and they were randomly selected from the test data. In the test data generation, standard deviation of model parameter distribution was assumed to be 20% of the standard values. The estimated hemoglobin changes are shown as light red (HbO) and light blue (HbR) solid lines. The assumed original hemoglobin changes are also shown in the figures.

064025_1_037906jbo6.jpg

Figure 6c is the estimation result where we applied an ordinary single-detector NIRS algorithm to the test data under the same conditions. The source-detector separation was set to be 30mm in this case. The sum of the partial path lengths of all layers was considered as the total path length in this case. The obtained hemoglobin changes were significantly smaller than the original changes because the total path length used in the algorithm is longer than the partial path length of the gray-matter layer. This is called partial volume effect.20

Next, we normalized the scale of each response in Figs. 6a, 6b, 6c by making its value at t=15s , 1 for HbO and 1 for HbR, to analyze the response shapes without being affected by their scale variation. The original changes were also normalized in this way. The results are shown in Figs. 6d, 6e, 6f. Figures 6d and 6e show that the shapes of the responses estimated by multidistance measurement are close to the original hemoglobin changes in the gray-matter layer, although the rising and decaying phases of the response are distorted. These distortions are considered to be caused by cross talk.15, 20 On the other hand, the shapes estimated by an ordinary single-detector NIRS algorithm [Fig. 6f] are different from the original changes in gray-matter layer but close to those in the superficial layer. These results show that the multidistance measurement is effective to extract the correct hemoglobin changes in the gray-matter layer.

5.

Conclusion

The multidistance measurement method can separate the functional response of brain activity from the systemic physiological interference caused by task execution. However, because measurement weighting is difficult to determine experimentally, we estimated it by a Monte Carlo simulation. A five-layered slab model was used to represent a human adult head. Because model parameters such as the thickness and transport scattering coefficient depend on a subject and a measurement position and are difficult to measure, we assumed that model parameters obey normal distributions around standard parameter values. The measurement weighting that produces a good average separation performance for the model parameter distribution was determined.

Partial path lengths of superficial and gray-matter layers were computed based on standard parameters by a Monte Carlo simulation, and their sensitivity coefficients, which indicate how much the partial path length changes when model parameter changes, were also given. The results indicate that the partial path length of the superficial layer is much longer than that of the gray-matter layer. The thickness parameter of scalp and skull layers has a significant influence on the partial path length of the gray-matter and superficial layers. The transport scattering coefficient of the skull layer has also a comparable influence on the partial path length of the gray-matter layer. Partial path lengths are insensitive to all other parameters.

Using the obtained partial path lengths and their sensitivity coefficients, a measurement weighting designed to yield the average optimum separation performance for model parameter distribution was derived. The obtained optimum weighting is robust to model parameter deviation and realizes smaller MSEs compared to a weighting that is determined without considering model parameter distribution.

When optimizing a measurement weighting, the covariance matrix Σx of hemoglobin changes is required. In this paper, we assumed a simple model of a hemodynamic change time course according to the recent study.18 Σx was calculated by this model. Although we considered a statistical distribution of layer parameters, the variations of hemodynamic change was not considered in this paper to simplify the analysis.

A significant systemic physiological interference, which is highly correlated with the functional response, is observed for a motor task (finger opposition), and only a small interference is observed for a cognitive and a visual tasks.5 Thus, the proposed method is effective mainly for a motor task experiment. However, this method can also suppress the systemic interference in the extracerebral layer caused by heart beat and respiration, although the theoretical analysis is given without considering these interferences in this paper to make the analysis simple. Thus, the proposed method is expected to be effective for many cases.

In our method, the hemoglobin concentration change is estimated by applying an inverse of molar absorption coefficient matrix and an optimized measurement weighting to the absorbance change obtained by the multidistance measurement [Eq. 13]. The applicability of the method strongly depends on the accuracy of the assumed model and its parameters. Thus, its naive application to a real measurement may be inappropriate. As shown in Fig. 6, however, because an ordinary single-detector NIRS measurement may lead to large estimation errors if there exists a hemoglobin concentration change in a superficial layer, we think that the proposed method is valuable even when we have only a rough knowledge of the model and its parameter values.

Appendices

Appendix: Derivation of Optimum Measurement Weighting, W

The criterion J of Eq. 20 is represented as follows:

Eq. 27

J=tr{FWTGW}2tr(HWT)+trΣx,
where

Eq. 28

F=1Kk=1K(LkEΣxETLkT),

Eq. 29

G=(E1)TE1,

Eq. 30

H=(E1)TΣxET(1Kk=1KLkT).
First we decompose F , G , and H into block matrices according to the decomposition structure of W .

Eq. 31

F=(F11F12F21F22),

Eq. 32

G=(G11G12G21G22),

Eq. 33

H=(H11H12H21H22).
Using these Eqs. 31, 32, 33, Eq. 27 becomes

Eq. 34

J=tr{F11W1TG11W1}+tr{F12W2TG21W1}+tr{F21W1TG12W2}+tr{F22W2TG22W2}2tr(H11W1T)2tr(H22W2T)+trΣn.
Then, we have

Eq. 35

JW1=2G11W1F11+G21TW2F12T+G12W2F212H11

Eq. 36

=2G11W1F11+2G12W2F212H11,
where we rely on the property that F , G , and Σn are symmetric. Similarly,

Eq. 37

JW2=2G21W1F12+2G22W2F222H11.
Thus, we have

Eq. 38

G11W1F11+G12W2F21=H11,

Eq. 39

G21W1F12+G22W2F22=H11.
We write csA as a vector obtained by aligning column vectors of a matrix A in a vertical order. Then the preceding equations become as follows using Kronecker products:

Eq. 40

{F11G11}csW1+{F21TG12}csW2=csH11,

Eq. 41

{F12TG21}csW1+{F22G22}csW2=csH11.
Summarizing equations 40, 41, we have

Eq. 42

({F11G11}{F21TG12}{F12TG21}{F22G22})(csW1csW2)=(csH11csH22).
Thus, the optimum W1 and W2 are obtained as follows:

Eq. 43

(csW1csW2)=({F11G11}{F21TG12}{F12TG21}{F22G22})1(csH11csH22).

Acknowledgment

The authors thank the Photon Migration Imaging Laboratory at Massachusetts General Hospital for the Monte Carlo simulation code (available at http://www.nmr.mgh.harvard.edu/PMI/index.htm).

References

1. 

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 044014 (2007). https://doi.org/10.1117/1.2754714 1083-3668 Google Scholar

2. 

Q. Zhang, E. N. Brown, and G. E. Strangman, “Adaptive filtering to reduce global interference in evoked brain activity detection: a human subject case study,” J. Biomed. Opt., 12 064009 (2007). https://doi.org/10.1117/1.2804706 1083-3668 Google Scholar

3. 

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 1874 –1882 (2005). https://doi.org/10.1364/JOSAA.22.001874 0740-3232 Google Scholar

4. 

M. A. Franceschini, S. Fantini, J. H. Thompson, J. P. Culver, and D. A. Boas, “Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging,” Psychophysiology, 40 548 –560 (2003). https://doi.org/10.1111/1469-8986.00057 0048-5772 Google Scholar

5. 

M. A. Franceschini, D. K. Joseph, T. J. Huppert, S. G. Diamond, and B. A. Boas, “Diffuse optical imaging of the whole head,” J. Biomed. Opt., 11 054007 (2006). https://doi.org/10.1117/1.2363365 1083-3668 Google Scholar

6. 

Y. H. Zhang, D. H. Brooks, M. A. Franceschini, and D. A. Boas, “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt., 10 011014 (2005). https://doi.org/10.1117/1.1852552 1083-3668 Google Scholar

7. 

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 062111 (2007). https://doi.org/10.1117/1.2814249 1083-3668 Google Scholar

8. 

M. Hiraoka, M. Firbank, M. Essenpresis, M. Cope, S. R. Arridge, P. van der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol., 38 1859 –1876 (1993). https://doi.org/10.1088/0031-9155/38/12/011 0031-9155 Google Scholar

9. 

Y. Fukui, Y. Ajichi, and E. Okada, “Monte Carlo prediction of near-infrared light propagation in realistic adult and neonatal head models,” Appl. Opt., 42 2881 –2887 (2003). https://doi.org/10.1364/AO.42.002881 0003-6935 Google Scholar

10. 

F. Fabbri, A. Sassaroli, M. E. Henry, and S. Fantini, “Optical measurements of absorption changes in two-layered diffusive media,” Phys. Med. Biol., 49 1183 –1201 (2004). https://doi.org/10.1088/0031-9155/49/7/007 0031-9155 Google Scholar

11. 

T. S. Leung, C. E. Elwell, and D. T. Delpy, “Estimation of cerebral oxy- and deoxyhemoglobin concentration changes in a layered adult head model using near-infrared spectroscopy and multivariate statistical analysis,” Phys. Med. Biol., 50 5783 –5798 (2005). https://doi.org/10.1088/0031-9155/50/24/002 0031-9155 Google Scholar

12. 

S. J. Matcher, C. E. Elwell, C. E. Cooper, M. Cope, and D. T. Delpy, “Performance comparison of several publishing tissue near-infrared spectroscopy algorithms,” Anal. Biochem., 227 54 –68 (1995). https://doi.org/10.1006/abio.1995.1252 0003-2697 Google Scholar

13. 

E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal,” Appl. Opt., 42 2915 –2922 (2003). https://doi.org/10.1364/AO.42.002915 0003-6935 Google Scholar

14. 

G. M. Hale and M. R. Querry, “Optical constants of water in the 200-nmto200-μm wavelength region,” Appl. Opt., 12 555 –563 (1973). https://doi.org/10.1364/AO.12.000555 0003-6935 Google Scholar

15. 

K. Uludag, M. Kohl, J. Steinbrink, H. Obrig, and A. Villringer, “Cross talk in the Lambert-Beer calculation for near-infrared wavelengths estimated by Monte Carlo simulation,” J. Biomed. Opt., 7 51 –59 (2002). https://doi.org/10.1117/1.1427048 1083-3668 Google Scholar

16. 

S. J. Matcher, M. Cope, and D. T. Delpy, “In vivo measurements of the wavelength dependence of tissue-scattering coefficients between 760 and 900nm measured with time-resolved spectroscopy,” Appl. Opt., 36 386 –396 (1997). https://doi.org/10.1364/AO.36.000386 0003-6935 Google Scholar

17. 

M. S. Cohen, “Real-time functional magnetic resonance imaging,” Methods, 25 201 –220 (2001). https://doi.org/10.1006/meth.2001.1235 1046-2023 Google Scholar

18. 

S. Boden, H. Obrig, C. Kohncke, H. Benav, S. P. Koch, and J. Steinbrink, “The oxygenation response to functional simulation: is there a physiological meaning to the lag between parameters?,” Neuroimage, 36 100 –107 (2007). https://doi.org/10.1016/j.neuroimage.2007.01.045 1053-8119 Google Scholar

19. 

T. J. Huppert, R. D. Hoge, S. G. Diamond, M. A. Francheschini, and D. A. Boas, “A temporal comparison of BOLD, ASL, and NIRS hemodynamic response to motor stimuli in adult humans,” Neuroimage, 29 368 –382 (2006). https://doi.org/10.1016/j.neuroimage.2005.08.065 1053-8119 Google Scholar

20. 

G. Strangman, M. A. Francheschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage, 18 865 –879 (2003). https://doi.org/10.1016/S1053-8119(03)00021-1 1053-8119 Google Scholar
©(2009) Society of Photo-Optical Instrumentation Engineers (SPIE)
Shinji Umeyama and Toru Yamada "Monte Carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy," Journal of Biomedical Optics 14(6), 064025 (1 November 2009). https://doi.org/10.1117/1.3275466
Published: 1 November 2009
JOURNAL ARTICLE
10 PAGES


SHARE
Advertisement
Advertisement
Back to Top