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

Functional Near-Infrared Spectroscopy (fNIRS) is a promising method to study functional organization of the prefrontal cortex. However, in order to realize the high potential of fNIRS, effective discrimination between physiological noise originating from forehead skin haemodynamic and cerebral signals is required. Main sources of physiological noise are global and local blood flow regulation processes on multiple time scales. The goal of the present study was to identify the main physiological noise contributions in fNIRS forehead signals and to develop a method for physiological de-noising of fNIRS data. To achieve this goal we combined concurrent time-domain fNIRS and peripheral physiology recordings with wavelet coherence analysis (WCA). Depth selectivity was achieved by analyzing moments of photon time-of-flight distributions provided by time-domain fNIRS. Simultaneously, mean arterial blood pressure (MAP), heart rate (HR), and skin blood flow (SBF) on the forehead were recorded. WCA was employed to quantify the impact of physiological processes on fNIRS signals separately for different time scales. We identified three main processes contributing to physiological noise in fNIRS signals on the forehead. The first process with the period of about 3 s is induced by respiration. The second process is highly correlated with time lagged MAP and HR fluctuations with a period of about 10 s often referred as Mayer waves. The third process is local regulation of the facial SBF time locked to the task-evoked fNIRS signals. All processes affect oxygenated haemoglobin concentration more strongly than that of deoxygenated haemoglobin. Based on these results we developed a set of physiological regressors, which were used for physiological de-noising of fNIRS signals. Our results demonstrate that proposed de-noising method can significantly improve the sensitivity of fNIRS to cerebral signals.


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., 2009Zhang et al., , 2012bTanaka 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 detrending 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).
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(Drummond, , 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., 2010Li et al., , 2013Zhang 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 timedomain 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 timefrequency representations of a signal with optimal resolution for Frontiers in Human Neuroscience www.frontiersin.org December 2013 | Volume 7 | Article 864 | 2 each frequency band. Continuous wavelet transform of a time dependent signal S(t) is defined as (Mallat, 2009): Where t and t 0 are time and time shift, respectively, α is a scale (of dimension time), and W S (α, t 0 ) 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 = t − t 0 α 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 0 ) is a complex function of scale α and position t 0 . 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 1 √ α 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 pseudofrequency simplifies to f a = 1 α . 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 0 ) 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 FIGURE 1 | (A) Scheme of the wavelet transform. The synthetic signal S 1 (t) (left) is convolved with complex mother wavelet function stretched with the scale parameter α (middle) to provide 2D wavelet scalogram on the right side. The real part of mother Morlet wavelet is plotted with red and imaginary part plotted with a blue line. Sign * indicates cross-correlation operation. The magnitude of the complex 2D wavelet scalogram is shown in color code. The synthetic curve S 2 (t) contains three harmonic components with the frequencies 0.3, 0.1, and 0.033 Hz, respectively, as well as additional white noise. The scalogram on the right hand side clearly demonstrates three components at scales 3, 10, and 34 s corresponding to the inverse frequencies of the signal components. (B) Example of a 2D wavelet coherence calculated for two synthetic signals with three coherent components. WCT S1 on the left and WCT * S2 in the middle are complex 2D wavelet scalograms of two signals (plotted as color coded magnitude values), both containing phase shifted coherent components with frequencies 0.3, 0.1, and 0.034 Hz, respectively, and additional non-coherent white noise. Sign × indicates scalar product operation. On the right side the magnitude value of wavelet coherence WCOH(S 1 , S 2 , α) is shown. Three peaks on the WCOH plot correspond to the three coherent components in the two signals in R-, M-, and A bands.

Frontiers in Human Neuroscience
www.frontiersin.org December 2013 | Volume 7 | Article 864 | 3 (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 1 (t) and S 2 (t) is defined as follows (Torrence and Webster, 1999;Grinsted et al., 2004): indicates the covariance or scalar product of the wavelet coefficients, W S 1 (α, t 0 ), and W S 2 (α, t 0 ) of signals S 1 and S 2 at scale α. R X, X (W S 1 , α) and R X, X (W S 2 , α) denote the autocorrelation or power of W S 1 (α, 0) and W S 2 (α, 0). In other words, the complex function WCOH(α) represents the crossspectral 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 1 (t) and S 2 (t) are hypothetical physiological signals containing coherent oscillations in R-, M-, and A-bands. From left to right 2D wavelet scalograms of S 1 (t) and S 2 (t) are plotted together with the wavelet coherence (outmost right plot). The first signal S 1 (t) is identical to the example provided in Figure 1A. The second signal S 2 (t) also contains white Gaussian noise added to the sinusoidal components. S 2 (t) contains the same frequencies as S 1 (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 1 (t) and S 2 (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(Wabnitz et al., , 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 0 , corresponding to the total photon count, (ii) the 1st moment of DTOF, m 1 , 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 0 are analogous to signals measured in conventional continuous wave (cw) fNIRS experiments. HbO and HbR based on m 1 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 emitterdetector 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.

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 Frontiers in Human Neuroscience www.frontiersin.org December 2013 | Volume 7 | Article 864 | 5 FIGURE 3 | Scheme of wavelet coherence analysis (see Section Wavelet Coherence Analysis for details). In the first step fNIRS and physiological signals are wavelet transformed to provide wavelet scalograms WCT HbO, HbR and WCT phys . The coherence between these scalograms yields a scale dependent wavelet coherence. Frequency bands corresponding to the maximum of the wavelet coherence and the coherence phase at maximum are determined. These parameters are further used to create auxiliary physiological regressors in the GLM de-noising (see Section General Linear Modeling and Physiological De-noising and Figure 4).
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 0 -based HbO and HbR signals. The noise amplitudes were chosen in a way, that the photon noise in the synthetic m 0 -based signal matched the level of photon noise in m 1 -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 1and 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. 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 In the first step physiological signals are filtered and time shifted using frequency bands and time shifts determined in the group wavelet coherence analysis as described in Section Section Wavelet Coherence Analysis and illustrated in Figure 3. In this way the set of auxiliary physiological regressors is created. This set is completed with predictors modeling cerebral activation and is used in the GLM analysis (see Section General Linear Modeling and Physiological De-noising and this figure).

WAVELET COHERENCE ANALYSIS
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 pseudofrequency 0.029 Hz in A-band) is only detectable for wavelet coherence between m 0 -based HbO and SBF.
Wavelet coherence values obtained between all physiological signals and HbO based on different moments are clearly different, with m 0 -based signals exhibiting the highest correlation at all scales, followed by m 1 -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 0 -and m 1 -based signals. Interestingly, the V-based HbO and HbR wavelet coherence values are similar at all scales, while the m 0 -and m 1 -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 0 -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 1 -and V-based signals may be due to higher photon noise levels as compared to m 0 -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 0 -based signal. The dashed lines represent the same signal, but numerically matched to the level of photon noise present in m 1 -and V-based signals. A slight influence of photon noise is clearly visible for low scales in the Rband. 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  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 1 -and V-based signals in comparison with the m 0 -based signal (see Table 2). The difference in time delay between m 0 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 0 -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 Mband. 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 Frontiers in Human Neuroscience www.frontiersin.org December 2013 | Volume 7 | Article 864 | 8

FIGURE 7 | Magnitude values of group average of wavelet coherence between haemoglobin concentration changes HbO (left) and HbR (right) and the four physiological signals.
In the columns from top to bottom results for mean arterial blood pressure, heart rate, skin blood flow, and respiration are shown, respectively. Haemoglobin concentration changes extracted from m 0 , m 1 , and V are plotted in blue, green, and red, respectively. The gray area under the black curve indicates the noise level for the magnitude of the group mean value of the WCOH between HbO or HbR and the physiological signals. Light blue, light, red and light green bars indicate the following three time scale ranges: R-band (scales around 3 s and pseudofrequencies around 0.3 Hz), M-band (scales around 10 s and pseudo-frequencies around 0.1 Hz) and A-band (scales around 35 s and pseudo-frequencies around 0.033 Hz), respectively. The colored shadowed areas represent standard error of mean for each curve. 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.

A-band
HbO, s HbR, s SBF m 0 (−6.9 ± 3.9) (*, T = −6,7) (0 ± 8) 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 0 -, m 1 -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 0 -based HbR is affected. The physiological noise related to the SBF changes is present only in m 0 and m 1based 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 0 -based HbO. Significant activation was identified for channel 4 only, for the m 1 -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.5 HZ)
The R-band correlation peak is most obvious for the wavelet coherence between m 0 -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 0 -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 0 -based signals. Significantly lower values were observed for m 1 -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 1 -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.15 Hz)
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  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 0 -, m 1 -, 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 0 , m 1 , 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 1 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 0 -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 denoising 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.04 Hz)
As one can see in Figure 7 the WCA reveals synchronous oscillations in the very low frequency band between the m 0 -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 0 -based HbO, and in the two medial channels in m 1 -based HbO and m 0 -based HbR and in one lateral channel in m 1based HbR. Since m 0 is much more sensitive to superficial tissue then m 1 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 0 -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 denoising 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 2 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 Frontiers in Human Neuroscience www.frontiersin.org December 2013 | Volume 7 | Article 864 | 14 method was shown to significantly improve the sensitivity of fNIRS to cerebral activation.