Identifying and quantifying main components of physiological noise in functional near infrared spectroscopy on the prefrontal cortex

Introduction

Functional Near-Infrared Spectroscopy (fNIRS) is a powerful tool to study functional organization of the prefrontal cortex (Scholkmann et al., 2013c). Due to absence of hair on the forehead, and relatively short distance between the forehead surface and the frontal cortex (Okamoto et al., 2004) this cortical area is well accessible by near-infrared light. However, despite these beneficial biophysical circumstances physiological noise generally limits overall fNIRS sensitivity and specificity on the prefrontal cortex (Tachtsidis et al., 2009; Aletti et al., 2012; Gagnon et al., 2012; Kirilina et al., 2012). Physiological fNIRS noise is induced by fluctuations in both blood volume and blood oxygenation in the extra- and intracranial tissues. Global and local blood flow regulation processes in these anatomical compartments lead to oscillations on multiple time scales. As a result fNIRS sensitivity to functional neuronal signals on the single subject level deteriorates and additional variance is added on the group level, due to inter-subject variability of the systemic and skin physiology.

In literature a number of methods were proposed for the separation of physiological noise from the cerebral activation (Saager and Berger, 2005; Katura et al., 2008; Gregg et al., 2010; Saager et al., 2011; Aletti et al., 2012). Comprehensive review of these methods can be found in Scholkmann et al. (2013c). Here we briefly summarize some of them. The majority of these approaches may be subdivided into three classes. Among them the most powerful methods are based on the idea of superficial signal regression. These methods require additional fNIRS channels with short spatial source-detector separation (Saager and Berger, 2005; Gregg et al., 2010; Saager et al., 2011; Gagnon et al., 2012; Funane et al., 2013). Nevertheless, these approaches still fail in removing physiological noise components originating from the cerebral compartment. Secondly, transient differences of functional signals and physiological noise have been exploited by a variety of methods (Kohno et al., 2007; Zhang et al., 2009, 2012b; Tanaka et al., 2013). These methods typically assume statistical independence between physiological noise and cerebral signal and fail to separate task-evoked responses of physiological parameters. In the third class, separation of forehead noise from cerebral fNIRS signals is achieved by exploring additional information from concurrent recordings of global systemic physiology, e.g., blood pressure and heart rate (HR) or local physiological parameters, e.g., skin blood flow (SBF) (Franceschini et al., 2006; Rowley et al., 2007; Tachtsidis et al., 2010; Minati et al., 2011; Patel et al., 2011; Takahashi et al., 2011; Sato et al., 2013). However, this separation is challenging, since the precise physiological fNIRS noise mechanisms are unknown. This shortcoming results from the lack of robust physiological models capable of relating physiological parameters like e.g., blood pressure or SBF with the fNIRS observable, tissue haemoglobin concentration. Further complications arise from the complex structure of physiological noise. Physiological noise involves components originating from both cerebral and extra-cerebral tissue, and global systemic as well as local tissue specific regulatory processes. Moreover different physiological processes dominate at different time scales.

The goal of the present study is to improve this situation by identifying the global systemic and local tissue-specific physiological noise processes in fNIRS recordings at the forehead and to quantify their relative impact in a group analysis.

Below we shortly review the literature on physiological processes which could potentially contribute to the physiological noise in fNIRS at different time scales. The most prominent, however, unproblematic, component is the heart pulsation with frequencies around 1 Hz, which may be effectively removed by low-pass filtering, since the corresponding frequency band is well distinguished from that of hemodynamic responses. Simple de-trending or high pass filtering also removes very slow drifts on the time scale of several minutes. The major concerns for fNIRS data quality are oscillations in the range from 0.2 to 0.005 Hz, since they account for the main components of physiological noise. Moreover, these fluctuations may in the worst case be synchronized with the task and might thereby lead to false positives in fNIRS activation maps (Tachtsidis et al., 2009).

Based on physiological considerations it is common to distinguish three distinct bands: high frequency oscillations (HFOs) ranging from 0.2 to 0.5 Hz, low frequency oscillations (LFOs) from 0.08 to 0.15 Hz and very low frequency fluctuations (VLFOs) from 0.02 to 0.08 Hz. In the following we will first focus on global systemic and then on local regulatory processes specific to cerebral or to extra-cerebral compartments.

HFOs in fNIRS signals, with frequencies around 0.3 Hz, are predominantly induced by direct or mediated influence of respiration. In the following we therefore refer to the frequency band around 0.3 Hz as R-band.

High variability in global systemic parameters of blood pressure and HR was found in the frequency band of LFOs (Julien, 2006). First reported by Mayer in 1876 (Mayer, 1876), these slow blood pressure waves are usually referred to as Mayer waves and the corresponding frequency band as M-band. In how far Mayer waves are propagated into cerebral and skin compartments remains controversial (Tong and Frederick, 2010; Tong et al., 2011).

In the cerebral compartment LFOs, independent from blood pressure fluctuations, were observed with Laser Doppler flowmetry (LDF) (Hudetz et al., 1992) and optical imaging in animals (Mayhew et al., 1996) and humans (Rayshubskiy et al., 2013) and with fNIRS on the human cortex (Elwell et al., 1999; Obrig et al., 2000; Tachtsidis et al., 2004).

Also SBF exhibits LFOs and VLFOs. LFOs in the M-band were observed even in isolated vessels (Johansson and Bohr, 1966) and were shown to be related to the activity of vessel walls. The SBF fluctuations in VLFO band originate most likely from sympathetic control of the peripheral vasculature (Kastrup et al., 1989; Söderström et al., 2003). The forehead skin vasculature constitutes a separate layer of sympathetic and parasympathetic innervations (Drummond, 1996, 1997). This fact allows for a synchronization between skin VLOFs and certain type of tasks in the forehead (Tachtsidis et al., 2009; Kirilina et al., 2012). Since the typical period of stimulation used in block design fNIRS studies often corresponds to VLFO frequency band, we refer to this band as A-band (activation) throughout this manuscript.

In summary, physiological fNIRS noise originates from different global systemic and local regulatory processes, each dominating at a different time scale. The goal of the present study was to identify the dominating processes and quantify their relative impact. We furthermore aim to develop a method to account for the contribution of each noise process in the fNIRS analysis, and therefore provide a method for physiological de-noising of fNIRS data. In order to achieve this challenging task we combined depth selective time-domain fNIRS measurements with concurrent recordings of global systemic peripheral physiological measurements of respiration, mean arterial blood pressure (MAP) and HR as well as a local SBF recordings on the forehead. Temporal disentanglement of the different physiological processes was realized by an advanced WCA. This method, originally developed in the field of geo science and meteorology (Torrence and Webster, 1999; Grinsted et al., 2004) but also applied to fNIRS signals analysis (Rowley et al., 2007; Li et al., 2010, 2013; Zhang et al., 2012a). WCA allows to target the coherent content of two temporal signals specifically at multiple time scales. Based on results of the WCA we developed a physiological de-noising method for fNIRS signals based on General Linear Modeling (GLM) (Kiebel and Holmes, 2007) and auxiliary physiological regressors (Tachtsidis et al., 2010). This approach allowed us to develop a novel versatile but yet robust physiological fNIRS de-noising method.

In the Theory section we briefly describe the theoretical background of wavelet transform and WCA used in this study to decompose physiological noise in fNIRS into components coherent with particular physiological processes. In the Methods section we describe the setup used for time-domain fNIRS experiments on the forehead. Based on the findings summarized in the Results section we discuss the possible origin of three processes contributing to the physiological noise and analyse the performance of the developed de-noising technique.

Theory

Wavelet Transform

Continuous wavelet transform allows for constructing time-frequency representations of a signal with optimal resolution for each frequency band. Continuous wavelet transform of a time dependent signal S(t) is defined as (Mallat, 2009):

Where t and t₀ are time and time shift, respectively, α is a scale (of dimension time), and W_S(α, t₀) represents the wavelet transform of signal S(t). ψ(t) is the mother wavelet function. In the present study a Morlet wavelet (Goupillaud et al., 1984) was used as mother wavelet function. This complex function may be expressed as product of a harmonic function and a Gaussian:

where f_c and f_b are dimensionless parameters, determining the wavelet center frequency and wavelet bandwidth, respectively, and $s=\frac{t-t_0}{\alpha }$ is a dimensionless variable. We used f_c = 1 and f_b = 3, a parameter set that provides optimal trade-off between time and frequency resolution. W_S(α, t₀) is a complex function of scale α and position t₀. If scale α is greater than one, the mother wavelet function ψ is stretched in the time domain. If α is smaller than one (but positive) the function is compressed in time. The scaling factor $\frac{1}{\sqrt{\alpha }}$ in Equation (1) is to ensure the energy normalization for all values of scale α.

From Equation (2) one can see that the relationship between wavelet scale α and pseudo-frequency f_a is given by:

In our case the relationship between the scale α and the pseudo-frequency simplifies to $f_a=\frac{1}{\alpha }$.

To exemplify the above described algorithm, the wavelet transform of a synthetic signal S(t) is outlined from left to right in Figure 1A.

The synthetic signal S(t) is a sum of three weighted sinusoidal signals with frequencies of 0.3 Hz (R-band), 0.1 Hz (M-band), and 0.034 Hz (A-band) and white Gaussian noise. The signal S(t) models a hypothetical physiological noise process with three components of different physiological origin.

The absolute values of W_S(α, t) are depicted as 2D wavelet scalogram in the right part of Figure 1A. This representation has clear maxima at the three scales 3, 10, and 33 s, corresponding to the inverse frequencies of the components present in the signals.

Please note that the width of the large-scale (low-frequency) component is proportionally larger, than that of the low-scale (high-frequency) component. This reflects an intrinsic wavelet transform property, which adjusts the frequency resolution to the frequency band.

Wavelet Coherence Analysis

WCA is a tool to evaluate the correlation of two time dependent signals separately at different time scales.

Wavelet coherence of two signals S₁(t) and S₂(t) is defined as follows (Torrence and Webster, 1999; Grinsted et al., 2004):

R_{X, Y}(W_S1, W_S2, α) indicates the covariance or scalar product of the wavelet coefficients, W_S1(α, t₀), and W_S2(α, t₀) of signals S₁ and S₂ at scale α. R_{X, X}(W_S1, α) and R_{X, X}(W_S2, α) denote the autocorrelation or power of W_S1(α, 0) and W_S2(α, 0). In other words, the complex function WCOH(α) represents the cross-spectral power in two time series as a fraction of the total power of the series. The absolute value of wavelet coherence WCOH(α) ranges from 0 to 1 with 1 indicating strongest correlation. The phase of the wavelet coherence provides information about the time lag between two signals. We define a time lag between two signals at scale α as the product of α and the wavelet coherence phase at this scale (expressed in radians).

Note that instead of conventionally performed smoothing in both time and scale dimensions (Torrence and Webster, 1999), we calculated only an average over the whole measurement period, but omitted smoothing in the scale dimension.

To illustrate the ability of WCOH(α) to detect coherent content in two signals at multiple time scales, Figure 1B provides an example of WCA of two synthetic signals with coherent components in three frequency bands. Two signals S₁(t) and S₂(t) are hypothetical physiological signals containing coherent oscillations in R-, M-, and A-bands. From left to right 2D wavelet scalograms of S₁(t) and S₂(t) are plotted together with the wavelet coherence (outmost right plot). The first signal S₁(t) is identical to the example provided in Figure 1A. The second signal S₂(t) also contains white Gaussian noise added to the sinusoidal components. S₂(t) contains the same frequencies as S₁(t) (0.3, 0.1, and 0.034 Hz), but with different phases and weights for the three components.

The wavelet coherence plot in Figure 1B shows three distinct maxima at 3, 10, and 34 s (equal to inverse frequencies of signal components), corresponding to three coherent components contained in the signals S₁(t) and S₂(t). In contrast to a simple correlation analysis WCOH(α) provides detailed frequency information about the coherent content in both signals.

In the current study we explore this property of WCA in order to analyze coherence between global systemic and local physiological processes and fNIRS signals at different time scales. In this way we are able to decompose physiological noise in fNIRS into components induced by several physiological processes.

Note the different widths of the maxima at the three different scales in the Figure 1B. Due to the wavelet transform properties, the frequency resolution degrades with the increasing scale α. However, the ratio between the resolution and α is independent of α.

Methods

fNIRS and concurrently recorded physiological signals were acquired as part of a comparative fNIRS/fMRI experiment previously reported in (Kirilina et al., 2012). Here we only briefly describe one part of experimental settings including fNIRS and peripheral physiological measurements, which were used in the analysis presented. A detailed description of subject population, stimulation paradigm, NIRS instrumentation and fNIRS data pre-processing can be found in (Kirilina et al., 2012).

Subjects

Fifteen healthy subjects (5 female, age 34.9 ± 7.2 years) took part in the present study. Due to technical reasons data from one subject was excluded from further analysis. During the experiment subjects were sitting in upright position in front of a computer monitor, while responding with right hand button presses. All subjects gave their informed consent to the experimental protocol, which was approved by a local ethics committee.

Stimulus

A variation of continuous performance task (CPT) combined with a semantic categorization task was used to achieve the cerebral activation in the frontal lobe (bilateral Brodmann Area 10). A series of German words representing either concrete or abstract categories were continuously presented to the subjects in semantic (sem-CPT) and control (word-CPT) tasks. In sem-CPT the subjects were instructed to press the left button if a concrete word was presented after an abstract word. In any other case, they should press the right button. In word-CPT the subjects were instructed to press the left button if one particular target word 1 (VORZUG, German for “preference”) followed target word 2 (KOFFER, German for “suitcase”), and right button in any other case. Both tasks were performed in nine 34.15 s long blocks in alternating order, interleaved by 31.11 s long baseline blocks. During the baseline blocks a fixation cross was presented in the middle of the screen. Pre- and post-baselines of 120 s were recorded for each subject.

Data Acquisition

Concentration changes in oxygenated and deoxygenated haemoglobin, in the following referred as ΔHbO and ΔHbR, respectively, were measured by the PTB time-domain optical brain imager (Wabnitz et al., 2005, 2010). This device provides three wavelengths 689, 797, and 828 nm. The laser power was split to obtain two sources. Diffusely reflected light was collected by four detection fiber bundles and detected by fast photomultipliers connected to a multi-board time-correlated single photon counting (TCSPC) system. Distributions of photon time of flight (DTOFs) were acquired with time bins of 24.4 ps width and at a 20 Hz rate. A set of one source fiber and two detection fiber bundles was placed on the left and right forehead, respectively, along the Fp1–Fp2 line defined by the international 10–20 system, as illustrated in Figure 2. A source-detector separation of 3 cm was chosen for all fNIRS channels. To exploit the potential of time-domain fNIRS for the separation of extracranial and cerebral signals, the measured DTOFs were analyzed in terms of statistical moments (Liebert et al., 2004). ΔHbO and ΔHbR time courses were extracted from changes in three different measures. These were (i) the 0th moment of DTOF, m₀, corresponding to the total photon count, (ii) the 1st moment of DTOF, m₁, i.e., the photon mean time of flight and (iii) the 2nd central moment of DTOF, the variance V. A detailed description of the procedure is given in Appendix 1 of (Kirilina et al., 2012). We would like to note that haemoglobin concentration changes based on m₀ are analogous to signals measured in conventional continuous wave (cw) fNIRS experiments. ΔHbO and ΔHbR based on m₁ and on the variance signal V have the advantage to be more sensitive to deeper and less sensitive to superficial absorption changes (Liebert et al., 2004; Wabnitz et al., 2005).

Peripheral Physiology

Alongside with fNIRS recordings, the following global and local physiological signals were recorded: HR (from electrocardiography ECG), blood pressure, scalp blood flow, and respiration chest movement.

Respiration (RSP) and ECG were recorded using a Nexus-10 system (Mind Media, The Netherlands). The respiration was recorded with a Hall-sensor based respiratory belt placed over the subject's lower ribs. ECG was recorded with a sampling rate of 256 Hz, using two electrodes placed on the upper right and lower left part of the chest, respectively, and a ground electrode placed on the upper left chest.

HR was defined for each heartbeat as the inverse time interval between two subsequent R-peaks in the ECG time trace. A PortaPress system (TNO TPD Biomedical Instrumentation) was used to continuously measure MAP at the left hand index finger. HR and MAP signals originally measured for time points corresponding to heartbeats were re-sampled onto the equidistant 1 Hz time grid.

Changes in SBF were recorded by a floLAB Laser Doppler Perfusion Monitor (Moor Instruments) operating at an emitter-detector distance of about 1 mm. The laser Doppler probe was placed on the right forehead 15 mm above the NIRS source (see Figure 2). Due to the close placement of fNIRS detectors and laser Doppler probe the light emitted by the floLAB device was detectable by the td fNIRS system. This light appeared as a constant uncorrelated background in the td-NIRS measurement and was subtracted during the data preprocessing procedure. The additional contribution of the background to photon noise was insignificant.

Data Analysis

Pre-processing

Pre-processing, wavelet analysis and GLM analysis of fNIRS and physiological data was performed by own software written in MATLAB (R2012a, Mathworks Inc.). All signals were first filtered by a low-pass filter, with a cutoff frequency of 0.8 Hz to remove the variability due to the cardiac cycle. In addition, a high-pass filter with a cutoff frequency of 0.008 Hz was used to remove very slow signal and baseline shift variations (Mitsis et al., 2004). In both cases fifth order Butterworth filters in forward and backward direction were used for further data filtering. After filtering, all recorded measurements were down sampled to 1 Hz.

Wavelet coherence analysis

After preprocessing of ΔHbO and ΔHbR signals from four channels and three moments as well as the time dependent physiological traces for each subject, WCA was employed to investigate the coherence between physiological noise in fNIRS and peripheral physiological traces.

We used a complex Morlet mother wavelet as defined by function cmor3-1, [see Equation (2)] in the MATLAB Wavelet Toolbox. A complex wavelet transform was performed on an equidistant scale grid ranging from 0.2 to 50 s with 0.2 s steps, corresponding to pseudo-frequencies from 5 Hz down to 0.02 Hz. For each subject, wavelet coherence was calculated between four physiological traces and twelve haemoglobin concentration changes measured for the four channels and based on three moments. Magnitude and phase of the group average was calculated in the next step for each pair of signals. Based on the observed maxima of the magnitude of the group wavelet coherence we identified three bands. Mean phase differences and corresponding time lags between each signal pair were calculated within each scale band. The above described procedure of WCA between fNIRS and physiological signals is illustrated as a flowchart in Figure 3.

Both fNIRS and physiological signals contain a certain amount of measurement noise. Therefore, the obtained wavelet coherence values include a stochastic component induced by experimental noise. The noise propagation in wavelet coherence is not linear due to the presence of the denominator in Equation (4). Therefore, proper statistical analysis of the group average wavelet coherence values and their comparison becomes a highly complex challenge. This situation is further complicated by the fact that fNIRS signals based on different moments have different photon noise contributions and that the physiological traces have different frequency spectra.

In order to overcome this non-linear problem, we numerically estimated the impact of different levels of photon noise and uncorrelated physiological noise components on the group wavelet coherence values. This was achieved by adding artificial noise to experimental data. In the following the numerical procedure to estimate the influence of these noise components is described in detail.

Influence of photon noise on wavelet coherence. In order to be able to correctly compare wavelet coherences obtained for fNIRS signals based on different moments, their different levels of photon noise were taken into account. In general, photon noise has increasing influence for haemoglobin concentrations based on higher order moments. We first theoretically estimated the standard deviation of the photon noise component in the measured ΔHbO and ΔHbR signals based on the three different moments. Theoretical and numerical details of this estimation are described in the Appendix. In a second step we performed numerical experiments, adding a controlled amount of synthetic white Gaussian noise to the (least noisy) measured m₀-based ΔHbO and ΔHbR signals. The noise amplitudes were chosen in a way, that the photon noise in the synthetic m₀-based signal matched the level of photon noise in m₁- and V-based ΔHbO and ΔHbR signals, respectively. The group wavelet coherence, between synthetic “noise matched” signals and physiological traces, was calculated and compared to the wavelet coherence obtained for experimentally measured m₁- and V-based signals.

Uncorrelated physiological noise and wavelet coherence. Furthermore, we numerically estimated the noise levels of the magnitude of group wavelet coherence obtained for uncorrelated physiological noise. For this purpose we calculated the wavelet coherence on the single subject level, between the fNIRS signals measured for each subject and physiological signals measured for different randomly chosen subjects. Additionally we circularly shifted the physiological traces with a random time shift, different for each subject in order to avoid inter-subject coherence effects potentially induced by the task. The wavelet coherences of individual subjects were subjected to a group analysis by calculating the absolute value of the group mean wavelet coherence value.

General linear modeling and physiological de-noising

As was described above, the proper statistical analysis of group wavelet coherence is mathematically challenging. Therefore, in this study we used the results of the group WCA only in a qualitative manner. Based on the maxima in wavelet coherence we identified the main processes contributing to physiological noise in fNIRS. We then used this information in order to create a set of auxiliary regressors for physiological noise modeling as described in the next session. Finally with the help of the GLM analysis informed by WCA we were able to statistically analyze the impact of each physiological process at the group level.

The GLM analysis was performed on the fNIRS signal time traces from the entire experimental sessions. Our model to analyze fNIRS signals included one regressor modeling task-related brain activation and auxiliary regressors modeling physiological noise. To create a regressor modeling task-related brain activation, boxcar functions of the task presentation were convolved with the standard haemodynamic response function as it is implemented in SPM8 software package (Friston and Stephan, 2007). The physiological noise regressors were constructed based on the results of the WCA.

For each physiological trace we determined a frequency band demonstrating maximal wavelet coherence with fNIRS signals. The time delay between fNIRS and physiological signals was determined in these bands, based on the phase of the group averages of wavelet coherence.

Physiological signals of the entire experimental session of each participant were then band-pass filtered at these frequency bands using a fifth order Butterworth filter in forward and backward direction. After filtering the signals were time shifted by the delay determined as described before. Filtered and time-shifted signals were used as regressors in the GLM analysis. The above described GLM de-noising procedure is illustrated as a flowchart in Figure 4.

Signal amplitudes corresponding to each regressor, obtained in the first level analysis for each subject, were subjected to a second level analysis, by calculating one-sample T-tests. For those values, which significantly differed from zero, with a significance level of p = 0.05 the group mean and standard deviation were calculated.

Results

Wavelet Coherence Analysis

Figure 5A shows time traces of the single channel ΔHbO and ΔHbR signals, as well as time traces of four physiological signals for one representative subject.

Figure 6 shows the 2D wavelet scalograms of ΔHbO and MAP and wavelet coherence between these two signals for one channel of a representative subject. In this particular case an increased coherence is observed for scales around 10 s, corresponding to the M-band.

Figure 7 depicts absolute values of the group average of the wavelet coherence between transient haemoglobin concentrations, e.g., (ΔHbO and ΔHbR) and the physiological signals MAP, HR, SBF, and RSP, respectively. The black solid curves in each subplot in Figure 7 indicate the level of coherency obtained between corresponding fNIRS signal and uncorrelated physiological noise as described in section Uncorrelated Physiological Noise and Wavelet Coherence. Therefore, all values within gray area under the black curve can be considered as not significant.

There are three distinct peaks of the wavelet coherence between fNIRS and physiological signals within three different frequency bands. The peak with the lowest scale around 3 s (corresponding to the pseudo-frequency of 0.33 Hz in R-band) is most clearly observed for the respiration signal, but also present in the wavelet coherence with MAP, HR, and SBF.

The peak with scale in the range of about 10 s (corresponding to the pseudo-frequency of 0.1 Hz in M-band) is strongly present in both global signals (MAP and HR) and is also detectable for SBF. In comparison to HR and SBF, MAP demonstrates higher correlation and a broader peak in the M-band, ranging up to scale values of 20–25 s.

The third and broadest peak at 34 s (corresponding pseudo-frequency 0.029 Hz in A-band) is only detectable for wavelet coherence between m₀-based ΔHbO and SBF.

Wavelet coherence values obtained between all physiological signals and ΔHbO based on different moments are clearly different, with m₀-based signals exhibiting the highest correlation at all scales, followed by m₁-based and then by V-based signals. In contrast, wavelet coherence values obtained for ΔHbR, based on all three moments, are not significantly different from each other. Generally higher wavelet coherence values were observed for ΔHbO in comparison to ΔHbR for m₀- and m₁- based signals. Interestingly, the V-based ΔHbO and ΔHbR wavelet coherence values are similar at all scales, while the m₀- and m₁-based ΔHbO coherence is always higher than that of ΔHbR for all signals (see Figure 7).

The lower wavelet coherence obtained for fNIRS signals based on higher moments may be rationalized by two different explanations. The higher coherence of m₀-based signals may originate from larger contributions of physiological fluctuations in the extra-cerebral skin tissue in these signals. On the other hand, the lower coherence in m₁- and V-based signals may be due to higher photon noise levels as compared to m₀-based signals. The estimation of latter effect is presented in Figure 8. It illustrates a hypothetical case assuming that the difference between signals based on different moments is only attributed to the different photon noise level. The blue solid lines show the experimentally obtained wavelet coherence for the m₀-based signal. The dashed lines represent the same signal, but numerically matched to the level of photon noise present in m₁- and V-based signals. A slight influence of photon noise is clearly visible for low scales in the R-band. In contrast, the impact of photon noise may be neglected for scales higher than 5 s (see Figure 8).

In the following we describe the time lags between physiological signals and fNIRS signals in the three main frequency bands exhibiting high impact of physiological noise. These results are presented for the R-, M-, and A-bands in Tables 1–3, respectively.

The time delay between the respiration signal and ΔHbO was around 1 s, with no significant difference between delays obtained for ΔHbO based on different moments. The time delays for ΔHbR based on the three moments were not significantly different from zero (see Table 1).

The delay between ΔHbO and MAP was around −0.6 s and was significantly shorter for the m₁- and V-based signals in comparison with the m₀-based signal (see Table 2). The difference in time delay between m₀ and V was (0.27 ± 0.1) s. The time delay between ΔHbR and MAP was around 3.3 s and was not significantly different for signals based on different moments.

The time delays between ΔHbO and SBF in the A-band are presented in Table 3. A significant delay of −6.0 s was observed for the m₀-based ΔHbO signal.

The results of group wavelet coherence between the four physiological signals are presented in Figure 9. A maximum in wavelet coherence in R-band can be seen for all pairs of physiological signals, with highest coherence in this band observed between MAP and RSP and between HR and RSP. Very high coherence (mean value 0.85) was observed between MAP and HR signals in the M-band. Maxima of lower amplitude in the M-band were observed for coherence between MAP and SBF as well as between HR and SBF. No considerable coherence was observed between any pair of physiological signals in the A-band. Table 4 shows the group average time delays between physiological signals in R- and M-bands.

Physiological De-Noising

Based on the qualitative analysis of group wavelet coherence presented in Figure 7, we identified the most relevant frequency bands and most relevant sources of physiological noise in each band. We then constructed auxiliary physiological regressors modeling impact of each process on fNIRS signals. Taking into account the mutual correlation of physiological signals that is obvious from the results presented in Figure 9, only a single physiological signal, the one showing the highest coherence, was used in each band for the de-noising procedure. In this way we avoid redundancy in our model. The RSP signals was used in order to model physiological noise in the R-band, the MAP signal was used in order to model the physiological noise component in the M-band. The physiological noise modeling in the A-band was based on the SBF signal.

In the following we summarize filter functions and time shifts employed for each physiological auxiliary regressor. This parameter choice was based on maxima in the group wavelet coherence (see Figure 7) and group mean time shifts for the corresponding bands (see Tables 1–3). To model respiration-induced noise, RSP signals were filtered with a bandpass filter [bandwidth (bw) 0.2 to 0.5 Hz] and time shifted by 1 and 0 s for ΔHbO and ΔHbR, respectively.

The impact of physiological noise at M-band frequencies was taken into account by bandpass filtering of MAP signals (bw 0.15–0.08 Hz) and a time shift of −0.69 s for ΔHbO and 3.4 s for ΔHbR signals.

To account for physiological noise in the A-band, SBF signals were band pass filtered (bw 0.02–0.04 Hz) and shifted in time by −7 s for both, ΔHbO and ΔHbR signals.

Time shifted and band-pass filtered physiological signals of each subject were normalized to unit power and used to create an individual set of four auxiliary physiological regressors for each subject. Figure 5B shows the time traces of four GLM regressors for one representative subject. GLM including one functional (cerebral) and four physiological regressors was performed for the four detector channels and two haemoglobin concentrations for each of the 14 subjects. In order to quantify the impact of physiological de-noising, we additionally performed GLM modeling with a reduced model, including the cerebral regressor only. The results of subject level analysis were subjected to a group T-test to determine the significance of activation at the group level.

The results of the second level analysis are presented in Table 5, and in Figures 10, 11 for ΔHbO and ΔHbR, respectively. In Figures 10A, 11A ΔHbO and ΔHbR attributed to cerebral activation are shown. Significant positive HbO and significant negative HbR concentration changes were observed in channel 4 for m₀-, m₁- as well as for V-based signals. For the V-based ΔHbR signal significant negative changes were also observed in channel 3.

Subpanels B–D in Figures 10, 11 show significant ΔHbO and ΔHbR values attributed to the three components of the physiological noise related to M-wave, SBF and respiration. One can see that the impact of physiological noise is generally much stronger for ΔHbO than for ΔHbR. Noise related to M-wave and respiration is present in ΔHbO based on all three moments, while essentially only m₀-based ΔHbR is affected. The physiological noise related to the SBF changes is present only in m₀ and m₁-based signals and is stronger pronounced in the two medial channels. In the V-based ΔHbR no significant contribution of physiological noise was detected.

The results of the reduced GLM including only one cerebral regressor are presented in Table 6. No significant cerebral activation was observed with this reduced model for m₀-based ΔHbO. Significant activation was identified for channel 4 only, for the m₁- and V-based ΔHbO signals. In ΔHbR signals significant activation was observed in channel 4 for the signals based on all three moments, however, the obtained T-values were always lower than that observed with the full model (compare Table 5, 6).

Discussion

The wavelet coherences of fNIRS signals and physiological processes, presented in Figures 7, 9, reveal a strong impact of the physiological noise at three main frequency bands: R-band with scales around 3 s (0.3 Hz), M-band with scales around 10 s (0.1 Hz) and A-band with scales around 34 s (0.034 Hz). These three distinct bands indicate the presence of at least three distinct physiological mechanisms dominating physiological fNIRS noise. In the following we summarize the results obtained at each band and discuss possible underlying physiological mechanisms.

R-Band (0.2 to 0.5Hz)

The R-band correlation peak is most obvious for the wavelet coherence between m₀-based ΔHbO signal and the respiratory signal (see Figures 7, 10D). Its peak frequency (0.3 Hz) corresponds to the mean respiration rate, thus indicating the direct influence of the respiration on the fNIRS signal. Clear coherence at the R-band was also found for m₀-based MAP, HR and SBF signals for both ΔHbO and ΔHbR concentration changes. However, the peak is not always clearly detectable, due to an overlay from neighboring M-band maxima in these signals.

Coherence at R-band demonstrates a clear difference between signals based on different moments. Highest coherence was primarily observed for m₀-based signals. Significantly lower values were observed for m₁-based ΔHbO and ΔHbR signals. The coherence in the R-band for V-based signal does not exceed the noise level. However, this dependence on the moment, which at first glance seems to be related to different physiological noise levels at different depths, has to be interpreted with care.

The analysis of the photon noise influence presented in Figure 8 indicates that the wavelet coherence at small scale values is sensitive to photon noise. Taking into account an increasing photon noise contribution in m₁- and V-based fNIRS signals, the observed dependence might be partly due to an artifact induced by different photon noise levels in the three measurements based on different moments.

In literature two different physiological mechanisms leading to R-band fluctuations are described. First, there is a direct influence of the respiration on the venous blood flow. Negative intrathoracic pressure during inspiration leads to increased venous outflow, which modulates the venous blood volume with the respiration frequency. The second mechanism is the respiratory modulation of the HR, mediated by the parasympathetic nervous system. Consequently the arterial input is modulated via HR by the respiratory cycle, a phenomenon usually referred to as respiratory sinus arrhythmia (Hirsch and Bishop, 1981). High correlation between respiration and HR observed in our study (see Figure 9) also indicates the presence of this phenomenon. Based on our data it is difficult to conclude, which of the two mechanisms dominates fNIRS physiological noise—direct influence of the respiratory pump, or the indirect influence of the HR variability. It is possible that the impact of both mechanisms is different for the two haemoglobin concentrations.

M-Band (0.05 to 0.15Hz)

The highest wavelet coherence between fNIRS signals and MAP and HR traces was observed in the M-band. This coherence is induced by Mayer waves, correlated fluctuation between HR and MAP in this frequency band (Elstad et al., 2011). As one can see in Figure 7, Mayer waves are present in both ΔHbO and ΔHbR concentration changes. However, both haemoglobin concentrations demonstrate different amplitudes, different time shifts and different dependence on moment, indicating different depth localization in the tissue. As one can see from the results of the WCA presented in Figure 7 and from GLM results in Figure 10, m₀-based ΔHbO signals show the highest M-band contribution. The amplitude of M-band physiological noise is higher than that of the cerebral signal. ΔHbO signals based on m₁ and V show significantly lower Mayer wave contributions. Since the influence of photon noise is negligible in the M-band, this moment dependence can be clearly assigned to Mayer waves with different amplitudes at different depths. A large fraction of the M-band contribution in the ΔHbO signal originates from superficial extra-cerebral tissue, however, there is an additional cerebral contribution as well. The presence of two separate M-band components (extra- and intracranial) is further supported by the phase difference detected between m₀-, m₁-, and V-based signals and reported in Table 2. Indeed, if the Mayer waves appear both in the brain and in the scalp, one might expect a time delay between these two compartments, caused by different vascular path lengths as well as possible delays in sympathetic mediating signals between the two compartments (Tong et al., 2011). Although cerebral auto-regulation mostly acts at lower frequencies (Latka et al., 2005; Rowley et al., 2007), there might be a time shift due to this vascular property of brain vessels as well. Since m₀-, m₁-, and V-based signals reflect a linear combination of signals from skin and cerebral compartments with different weights, they would then show different time lags relative to MAP.

Interestingly, wavelet coherence between ΔHbR and MAP and between ΔHbR and HR is very similar for m₀, m₁, and V-based signals. In addition, no significant time lag was obtained for these signals. This means that the main source of M-band physiological noise in ΔHbR is localized in deeper intracranial tissue. However, GLM analysis failed to detect significant ΔHbR M-band contribution in V-based signals (see Figure 11), although significant M-band signals were obtained in m₁ and V-based ΔHbO signals. Therefore, we can conclude that both extra- and intracranial M-band signals are mostly present on arterial side. This fact might explain why Mayer waves are not considered as important source of physiological noise in fMRI (Birn et al., 2006; Chang and Glover, 2009). Since BOLD signal exploited in fMRI mostly reflects changes in HbR concentration, it seems to be not strongly affected by Mayer waves. Moreover, the cerebral compartment might exhibit an additional contribution in the M-wave band which is not coupled with blood pressure fluctuations as demonstrated recently by (Rayshubskiy et al., 2013). The presence of an additional component that is uncorrelated with MAP could also explain the experimentally observed lower wavelet coherence in the deeper tissue.

Another interesting fact is the high coherence between fNIRS signals and HR time traces in the M-band. Although the direct correlation between HR and signals is low, it accounts for up to 40% fNIRS signal variance for conventional cw (in our case m₀-based) signals, for HR being shifted in time. Despite this number is lower than that obtained for MAP, it is still enough to significantly improve fNIRS sensitivity when used in a de-noising procedure. This may become a fact of high importance in fNIRS experimental practice. Indeed, continuous blood pressure measurements are challenging and dedicated hardware more rare and cost intensive than conventional ECG or pulse plethysmography devices. Those are relatively cheap and available in many labs. Thus, physiological de-noising based on HR measurements is much more easily applicable than approaches based on continuous MAP recordings, although it might be slightly less efficient.

High correlation between MAP, HR, and fNIRS signals was reported and emphasized by several studies (Franceschini et al., 2006; Katura et al., 2006; Minati et al., 2011; Li et al., 2013). The time shift between fNIRS signals and MAP and HbO signals reported by (Katura et al., 2006) are very close to those presented in the Table 2.

Finally we would like to discuss the possible physiological mechanism inducing Mayer waves in fNIRS signals. Synchronized oscillations with frequencies around 0.1 Hz are typically observed in blood pressure and HR in humans and animals and are most likely induced by an interplay between sympathetically driven HR variations and sympathetic vasoconstriction of peripheral resistive vessels (Pagani et al., 1986; Malliani et al., 1991; Stauss et al., 1998; Cohen and Taylor, 2002; Nilsson and Aalkjaer, 2003). The high coherence between MAP and HR in M-band observed in our experiment (see Figure 9) is in agreement with these findings.

A mechanism that directly links local blood pressure and local HbO and HbR concentration changes is vasoconstriction of peripheral resistance vessels. These vessels are situated prior to the capillary bed on the arterial side and therefore reveal high concentration of oxygenated haemoglobin. The diameter of resistive vessels is sympathetically regulated, and critically influences the overall hydrodynamic resistance of the vascular system and thereby the blood pressure. Regulatory changes in the vessel diameter are connected to changes of the arterial blood volume, which is the most probable source of the observed tissue haemoglobin concentration changes. These changes occur on the arterial side are further propagated to the venous side due to induced fluctuations in blood flow (Tong et al., 2011).

A-Band (0.02 to 0.04Hz)

As one can see in Figure 7 the WCA reveals synchronous oscillations in the very low frequency band between the m₀-based ΔHbO signal and SBF.

With the help of GLM analysis presented in Figures 10, 11 we detect a significant contribution of this signal in all four channels in m₀-based ΔHbO, and in the two medial channels in m₁-based ΔHbO and m₀-based ΔHbR and in one lateral channel in m₁-based ΔHbR. Since m₀ is much more sensitive to superficial tissue then m₁ and V, we can conclude that activation in the A-band is predominantly localized in the skin compartment, and is due to the local SBF regulation mechanisms.

The stronger influence of the SBF related artifact on the medial forehead is supported by our earlier findings from a comparative fMRI/fNIRS study (Kirilina et al., 2012). In this study we observed a task-evoked response in the two medial veins draining the forehead.

The scale of maximum coherence corresponds to the period of stimulation (34 s) used in our study. Therefore, we hypothesize, that observed SBF responses might be induced by task-induced cognitive stress.

VLFOs with frequencies from 0.021 to 0.052 Hz were previously reported in LDF and NIRS measurements of the skin (Li et al., 2010) and were linked to sympathetic control of the peripheral vasculature (Kastrup et al., 1989; Söderström et al., 2003). Since increased level of sympathetic activity might be induced by cognitive or emotional tasks used in fNIRS experiments, SBF oscillation in this band might be synchronized with the task. Therefore, special caution has to be taken in fNIRS signal analysis in order to separate these superficial skin signals from cerebral activation.

Physiological De-Noising of fNIRS Data

The results of GLM-modeling presented in Figures 10A, 11A demonstrate significant activation in the lateral left channel. An increase in HbO and decrease in HbR concentration was observed in this area when physiological de-noising was applied. This result is consistent with the results of an fMRI study performed using the same task and the same subject group (Kirilina et al., 2012). In this study activation in bilateral Brodman Area 10 was observed, lateralized on the right side (Talairach coordinates [32 48 11]). The left BA10 activation was too deep from the cortical surface to be detected by fNIRS.

Although ΔHbO demonstrates the signals of higher absolute value than that of ΔHbR, it is also much more strongly affected by physiological noise (compare Figures 10B–D, 11B–D). Moreover, for the m₀-based signal we failed to detect cerebral activation in ΔHbO when the reduced GLM model (Table 6) was used. This means that we would not detect any cerebral activation with conventional cw fNIRS, if not corrected for physiological noise. In contrast, ΔHbR signals show robust activation in all three moments even without physiological noise correction (see Figure 11).

The above described results demonstrate that ΔHbR signals are more reliable in detecting cerebral activation, and that the de-noising procedure developed in the current study can significantly improve the sensitivity of fNIRS to cerebral activation for ΔHbO.

Limitations of the Current Study and Outlook

One limitation of the present approach is that not all possible contributions to fNIRS physiological noise might be detected by our approach, but only those, which also manifest themselves in MAP, HR, and SBF. In particular, Lased Doppler Flowmetry measures only capillary SBF and not that of the large vessels. We also did not account for the possible impact of fluctuations in arterial CO₂ concentration (Scholkmann et al., 2013a,b). Thus, some important sources of signal variance might be missing in our consideration.

Moreover the wavelet analysis is a linear transformation method, thus with the present approach we might miss more complex interrelations between physiological parameters and fNIRS signals, which are not captured by a simple linear relationship. In our WCA we used single physiological trace at a time. In the future the multi-variable coherence and correlation analysis such as the canonical correlation analysis could provide more complete picture (Caicedo et al., 2013).

Another limitation of the proposed de-noising strategy is the necessity to use additional physiological sensors and monitors. In particular, continuous monitoring of the blood pressure as well as laser flowmetry are not widely available in most labs. However, we have shown that even applying widely available and cost efficient sensors for respiration and HR might significantly improve fNIRS data quality.

An important interesting step in the further investigation might be to compare and combine the presented method of physiological de-noising with methods based on anatomical localization of cerebral and extra-cerebral tissue such as superficial signal regression.

Conclusion

In the current study we investigated the impact of global systemic and local regulatory physiological processes on physiological noise in fNIRS measurements and developed a method for physiological de-noising of fNIRS data. Global systemic processes were quantified by measuring MAP, HR, and respiration. Local regulatory processes in the skin were measured by means of SBF recordings. WCA was employed to characterize the contribution of these physiological parameters on fNIRS noise, at different time scales. Time-domain measurements in combination with signals, based on different moments of DTOF, enabled us to obtain information on the depth localization of different physiological noise sources.

With the help of these methods we were able to identify three main mechanisms contributing to physiological noise in fNIRS signals at three time scales. Two of these processes were induced by global systemic physiology: Mayer waves and respiration. The third slow process was induced by local blood flow changes in the skin tissue. We show that HbO signals are more strongly affected by global processes in both, extra- and intra-cerebral compartments, and local SBF regulation, while HbR signals are less contaminated by extra-cerebral processes.

By means of GLM analysis and auxiliary physiological regressors we quantified the relative impact of each process. Moreover, we propose a de-noising algorithm and demonstrate its performance on a functional experiment on the forehead. The proposed method was shown to significantly improve the sensitivity of fNIRS to cerebral activation.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The research leading to these results has received funding from the European Community's Seventh Framework Programme [FP7/2007-2013] under grant agreement n FP7-HEALTH-F5-2008-201076, and funding from German Research Organization (DFG) in the framework of Cluster of Excellence “Languages of Emotion” EXC302. IT and NY were supported by the Wellcome Trust (088429/Z/09/Z). The authors are thankful to Dr. Angela Heine for development of the task used in this study.

References

Appendix

The concentration changes of HbO and HbR were determined based on three moments of DTOFs (m₀, m₁, and V) separately. The time traces of ΔHbO and ΔHbR related to the different moments exhibit different levels of non-physiological noise, due to the different influence of instrumental noise. In order to estimate the amplitude of this non-physiological noise in ΔHbO and ΔHbR time traces, in the first step we estimated the instrumental and photon noise induced standard deviation of the corresponding moments. In the second step, we investigated how this uncertainty is propagated through the calculation of ΔHbO and ΔHbR.

Step 1. It is assumed that the uncertainty of measurement is dominated by photon noise following a Poisson distribution. For this case the variances of the changes in moments are determined by (Liebert et al., 2003, 2012):

where ΔM₀ = −ln(1 − Δm₀/m₀) is the attenuation change, ΔM₁ = Δm₁, ΔM₂ = ΔV. M_{2, raw} = V_raw and M_{4, raw} are the variance (second central moment) and the fourth central moment of the DTOF, respectively, as directly obtained in the experiment, without eliminating the influence of the instrumental response function. m₀, M_{2, raw} and M_{4, raw} were estimated from the experimentally measured DTOFs for each subject and each wavelength. Based on the group-averaged values of these parameters we estimated the standard deviations of each moment according to Equations (A1)–(A3) (see Table A1).

Step 2. In order to obtain the concentration changes of oxy- and deoxyhaemoglobin, ΔHbO_n and ΔHbR_n, based on the n-th moment we need to solve the over-determined set of linear equations:

where ΔM_n(λ_i) represents the changes in the n-th moment (n = 0, 1, and 2) (definitions see above) for the i-th wavelength. S_n(λ_i) is a moment-specific sensitivity factor, ϵ_HbO(λ_i) and ϵ_HbR(λ_i) are the molar absorption coefficients of HbO and HbR, The sensitivity factors (for the homogeneous case) can be derived from moments of the measured DTOF by (Arridge et al., 1992; Steinbrink, 2000; Liebert et al., 2004):

where M₁ = m₁, M₂ = V, M₃—third central moment of the DTOF and c_m—speed of light in the medium. Note that S₀ is related to the differential pathlength factor (D_PF) by S₀ = D_PFr_sd with r_sd being the source-detector separation.

Equation (A4) can be rewritten as $\overrightarrow{\Delta M_n}=A_n\overrightarrow{\Delta C_n}$ where $\Delta C_n=\left(\begin{array}{l}\Delta {\text{HbO}}_n\\ \Delta {\text{HbR}}_n\end{array}\right)$ is the vector of concentration changes, $\overrightarrow{\Delta M_n}$ the vector of changes of the n−th moment for the three wavelengths and A_n the coefficient matrix on the right-hand side of Equation (A4). The solution of this system of equations in a least-squares sense is given by:

The variances (i.e., squares of standard deviations) of the concentration changes ΔHbO_n and ΔHbR_n are the diagonal elements of the covariance matrix:

Herein Z_n is the covariance matrix of the moment changes $\overrightarrow{\Delta M_n}$. Since photon noise is uncorrelated for the measurements at different wavelengths, Z_n is a diagonal matrix:

with σ (ΔM_n(λ_i)) being the standard deviation of the n-th moment at wavelength λ_i.

For the estimation of the standard deviations of the concentration changes we used the group-averaged values of the sensitivity factors together with the group-averaged standard deviations of moments obtained in Step 1, both listed in Table A1.