NICA: A Novel Toolbox for Near-Infrared Spectroscopy Calculations and Analyses

Introduction

Functional near-infrared spectroscopy (fNIRS) is a non-invasive optical method for detecting functional hemodynamic activity of the cerebral cortex. This method measures the regional changes of the oxygen concentration in the blood, more precisely the concentration changes of oxygenated (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb; Jöbsis, 1977; Wolf et al., 2007; Scholkmann et al., 2014). These concentration changes are coupled with neuronal activity (i.e., neurovascular coupling) and can therefore be associated with activation of the cortex in the investigated area during certain cognitive and motor tasks (Villringer et al., 1993; Leithner and Royl, 2014).

In order to measure the concentration changes of oxy-Hb and deoxy-Hb, near-infrared light is sent from and detected by optical sensors (sources and detectors), so called optodes. The near-infrared light can penetrate the cranial bone, and the photons are absorbed and scattered by different interactions with the cerebral tissue. The scattered photons follow a random path through the tissue, where parts of these photons get absorbed (e.g., by chromophores such as oxy-Hb and deoxy-Hb). Some of the photons are back-scattered and can be measured by detectors which are placed at a distance of a few centimeters (∼3 cm) around the source (Okada et al., 1997; Boas et al., 2004; Ferrari and Quaresima, 2012). Using the (modified) Beer–Lambert law, it is possible to calculate the concentration changes in oxy-Hb and deoxy-Hb according to the attenuation of the intensity of the near-infrared light between the source and the detector (Coyle et al., 2007). With these concentration changes, conclusions about metabolic changes in the investigated area can be drawn, which are associated with brain activity (Malonek and Grinvald, 1996; Wolf et al., 2002).

Aside from brain activity, the recorded signals of the hemodynamic responses may contain influences from physiological signals which do not originate from the cerebral cortex (systemic influences) (Coyle et al., 2004; Takahashi et al., 2011; Tachtsidis and Scholkmann, 2016). These signals include the pulsation of the heart, the respiratory patterns and low frequency oscillations from the blood pressure (BP; De Boer et al., 1986; Koepchen, 1991). The frequency of the pulse waves (heart rate, HR) typically lies between 1 and 2 Hz, the respiratory frequency (RF) is in a range from 0.2 to 0.4 Hz (Bauernfeind et al., 2014). The third order BP waves, also known as Mayer–Traube–Hering (MTH) waves, occur between 0.07 and 0.13 Hz (Pfurtscheller et al., 2010b). Since the time constant of hemodynamic activation responses is about 5–10 s (Elwell et al., 1999), the systemic influences, especially the MTH waves and the RF, may influence and superimpose the recorded cortical activation (Obrig et al., 2000; Bauernfeind et al., 2008).

A number of methods have been proposed to reduce the interfering physiological signals, such as different types of low-pass filters (Bauernfeind et al., 2008; Power et al., 2011), moving average filters (Soraghan et al., 2008), as well as pulse regression (Gratton and Corballis, 1995) and spatial filtering (Pfurtscheller et al., 2010a). Although these methods can be very useful to remove physiological noise, they can also cause a reduction of the hemodynamic responses itself (Bauernfeind et al., 2014). Therefore, it is required that the systemic influences in the frequency range of the fNIRS signal are detected and removed. However, most studies show that these influences and their removal are neglected (Leff et al., 2011). Bauernfeind et al. (2014) proposed three different signal processing methods that can help to reduce the systemic influences of the recorded fNIRS-signal: a common average reference (CAR) method, independent component analysis (ICA) and the use of transfer function (TF) models. Therefore, it is necessary to simultaneously record the physiological signals (BP, RF, and HR).

To perform a thorough analysis of fNIRS data, it is necessary to pre-process the data, remove physiological and technical artifacts, calculate the concentration changes of oxy-Hb and deoxy-Hb and generate output files and figures to visualize the recorded signals. There are freely available toolboxes or compact software packages to analyze fNIRS data, such as Homer2 (Huppert et al., 2009), NIRS-SPM (Watanabe et al., 2018), or the NIRS Brain AnalyzIR Toolbox (Santosa et al., 2018). However, these toolboxes and software packages do not offer artifact correction methods for fNIRS data using simultaneously recorded physiological signals. For this purpose, the MATLAB toolbox NICA was developed – a graphical user interface (GUI), which offers a user-friendly handling of processing and visualizing fNIRS data, including all the necessary signal processing steps and the correction of physiological and technical artifacts.

With NICA it is possible to load and pre-process the fNIRS data, apply correction methods for physiological artifacts (systemic influences), exclude noisy channels and/or trials, and set different options for visualizing the results of the analysis. This article focuses on the implemented signal processing and physiological artifact correction, the implementation and handling of the GUI, as well as on the requirements for the measurement system and data. Furthermore, an example analysis will be provided, showing the importance of proper artifact correction.

Materials and Methods

FNIRS Measurement Device and File Requirements

NICA was designed and optimized to analyze measurement data recorded with the NIRScout 1624 device and its corresponding recording software NIRStar (NIRx Medizintechnik GmbH, Berlin, Germany). This measurement system uses a multichannel continuous wave (CW) technique with a sampling rate of 3.91 Hz. It offers 16 sources, where near-infrared light is emitted at a wavelength of 760 and 850 nm, and 24 detectors. In order to perform an analysis, two file formats are required for the software: a header-file (HDR-file) and an extensible data format file (XDF-file).

The HDR-file provides basic information about the measurement setup (filename, date, number of sources and detectors, etc.) and is generated with the NIRStar recording software (NIRx Medizintechnik GmbH, Berlin, Germany), which is used to measure the fNIRS signals. This software is generally used and delivered with the NIRScout device.

The XDF-file includes all the recorded signals of a measurement: fNIRS data, electrocardiogram (ECG), HR, BP, RF as well as the paradigm markers. The signals are streamed and synchronized in time via lab streaming layer (LSL) applications (Kothe, 2014) and recorded with the LabRecorder (default recording program for the LSL). The fNIRS signals are recorded with two different wavelengths, 760 and 850 nm, using the CW approach. ECG and RF are measured with surface electrodes and a stretching sensor in a chest strap (Respiratory Effort Sensor, Pro-Tech Services, United States), respectively, at a sampling rate of 256 Hz, and connected via a biosignal amplifier (g.USBamp, Guger Technologies, Austria) to the LSL application. HR and BP are recorded with a finger and an upper arm cuff with the CNAP Monitor 500 (CNSystems Medizintechnik AG, Austria). The CNAP device measures the blood pressure continuously, noninvasively and beat-to-beat on the finger artery with a sampling rate of 200 Hz. The upper arm measurement is used to calibrate the finger cuff to deliver absolute blood pressure values. The blood pressure values are displayed at the integrated monitor and sent via custom-made sync-box to the computer. The paradigm markers are sent during the paradigm presentation and are essential for the timing of the measurement (to distinguish between trials and different conditions). The physiological signals (ECG, HR, BP, and RF) are required for the artifact correction of the systemic influences. The ECG and the BP can both be used to correct the MTH-waves. If the CNAP Monitor is used, there is no need to measure the ECG as well. Even though an analysis can be performed without physiological signals, it is not recommended, since a thorough correction of the physiological artifacts is not possible.

Signal Processing and Artifact Correction

The signal processing steps in NICA from the raw signals to the clean oxy-Hb and deoxy-Hb signals, as well as the physiological artifact removal methods, were introduced by Bauernfeind et al. (2014). This article focuses on the implementation and utilization of NICA; therefore, this section only roughly describes the usage and the parameters of the signal processing and artifact correction methods.

Removing 50 Hz Power Line Interference

In the first signal processing step, a filter is applied on the ECG and the respiratory signal to remove the 50 Hz power line interference. The used filter is an infinite impulse response (IIR) notch filter (band-stop) with a cut-off frequency of 50 Hz and a bandwidth of 1.43 Hz. The removing of the power line interference of the other signals (fNIRS, BP, and HR) is included in their measurement device.

Calculating Concentration Changes of oxy-Hb and deoxy-Hb

The NIRScout measurement device uses the CW approach, which detects the attenuation of light as their only measurement parameter. The attenuation (A) can be calculated with Eq 1 (Delpy and Cope, 1997):

$$A={\log }_{10}\left(\frac{I_0}{I}\right)=\left({\mathrm{\alpha }}_{Hb}\cdot c_{Hb}+{\mathrm{\alpha }}_{Hb{O}_2}\cdot c_{Hb{O}_2}\right)\cdot x\cdot l+K(1)$$

A is the attenuation of light (extinction), I₀ the intensity emitted from the source, I the intensity measured from the detector, α_HbO2 and α_Hb are the molar coefficients of extinction of oxy-Hb and deoxy-Hb, c_HbO2 and c_Hb the corresponding concentrations of the absorber (oxy-Hb and deoxy-Hb), x is the differential path length (DPF), l the geometric distance between source and detector, and K describes the tissue loss due to scattering. The DPF cannot be determined with the CW approach and has to be based on measured values. Hence, it is not possible to calculate absolute values of the concentration of oxy-Hb and deoxy-Hb, but rather concentration changes (in mM × mm) (Duncan et al., 1995).

The tissue loss can be seen as constant in a short time interval and therefore be neglected after differentiation:

$$\mathrm{\Delta }A={\log }_{10}\left(\frac{I\left(t_1\right)}{I\left(t_2\right)}\right)=\left({\mathrm{\alpha }}_{Hb}\cdot \mathrm{\Delta }{c}_{Hb}+{\mathrm{\alpha }}_{Hb{O}_2}\cdot \mathrm{\Delta }{c}_{Hb{O}_2}\right)\cdot x\cdot l(2)$$

Rearranging Eq 2 leads to the following expression:

$$\frac{\mathrm{\Delta }A}{x\cdot l}={\mathrm{\alpha }}_{Hb}\cdot \mathrm{\Delta }{c}_{Hb}+{\mathrm{\alpha }}_{Hb{O}_2}\cdot \mathrm{\Delta }{c}_{Hb{O}_2}(3)$$

The calculation of the concentration changes of oxy-Hb and deoxy-Hb requires the two wavelengths with which the fNIRS signal was recorded (760 and 850 nm). Eq 3 can be used for both wavelengths and combined by using matrix notations:

$$\left[\begin{array}{c}\hfill \frac{\mathrm{\Delta }{A}_{760nm}}{x_{760nm}\cdot l}\hfill \\ \hfill \frac{\mathrm{\Delta }{A}_{850nm}}{x_{850nm}\cdot l}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathrm{\alpha }}_{760nm,Hb}\mathit{\hspace{1em}}\hfill & \hfill {\mathrm{\alpha }}_{760nm,Hb{O}_2}\hfill \\ \hfill {\mathrm{\alpha }}_{850nm,Hb}\mathit{\hspace{1em}}\hfill & \hfill {\mathrm{\alpha }}_{850nm,Hb{O}_2}\hfill \end{array}\right]\cdot \left[\begin{array}{c}\hfill \mathrm{\Delta }{c}_{Hb}\hfill \\ \hfill \mathrm{\Delta }{c}_{Hb{O}_2}\hfill \end{array}\right](4)$$

After rearranging Eq 4, the concentration changes of oxy-Hb Δc_HbO2 and deoxy-Hb Δc_Hb can be calculated with Eq 5:

$$\left[\begin{array}{c}\hfill \mathrm{\Delta }{c}_{Hb}\hfill \\ \hfill \mathrm{\Delta }{c}_{Hb{O}_2}\hfill \end{array}\right]={\left[\begin{array}{cc}\hfill {\mathrm{\alpha }}_{760nm,Hb}\mathit{\hspace{1em}}\hfill & \hfill {\mathrm{\alpha }}_{760nm,Hb{O}_2}\hfill \\ \hfill {\mathrm{\alpha }}_{850nm,Hb}\mathit{\hspace{1em}}\hfill & \hfill {\mathrm{\alpha }}_{850nm,Hb{O}_2}\hfill \end{array}\right]}^{-1}\cdot \left[\begin{array}{c}\hfill \frac{\mathrm{\Delta }{A}_{760nm}}{x_{760nm}\cdot l}\hfill \\ \hfill \frac{\mathrm{\Delta }{A}_{850nm}}{x_{850nm}\cdot l}\hfill \end{array}\right](5)$$

Removing Baseline Drifts From oxy-Hb and deoxy-Hb Signals

A 6th order Butterworth high-pass filter is applied to remove the baseline drifts from the oxy-Hb and deoxy-Hb signals. The cut-off frequency is 0.005 Hz, with a passband ripple of 1 dB and a stopband attenuation of 30 dB. The baseline drifts can be described as low frequency noise that shifts the baseline of the signal, which makes it difficult to detect characteristic features (Kaur and Singh, 2011).

Physiological Artifact Correction

As in Bauernfeind et al. (2014), three different signal processing methods are used by the toolbox that can help to reduce the physiological artifacts of the recorded fNIRS-signal: CAR, ICA, and TF models.

Common average reference

The CAR method is a simple technique to reduce noise (of technical or physiological origin) in the recorded signal, which has already been used for the analysis of EEG data (Zhang et al., 2007; Matthews et al., 2008; Ludwig et al., 2009). For analyzing fNIRS data, the CAR method was first introduced by Pfurtscheller et al. (2010a). This method uses the fact that systemic influences interfere with the signal in all channels, and can therefore be reduced by using the mean of all channels:

$$N_i\left[n\right]=X_i\left[n\right]-\frac{1}{M}\sum _{j=1}^M{X}_j\left[n\right](6)$$

The new signal N_i[n], as described in Eq 6, is calculated by subtracting the mean of all channels X_j, j = 1 … M from each single channel X_i, i = 1 … M for each time point n.

Independent component analysis

Independent component analysis is based on the mixing of a linear transformation Matrix A and the independent sources S, which results in the recorded and artifact influenced fNIRS signals X (Eq 7):

$$X=AS(7)$$

The estimates (independent components, ICs) of the sources (oxy-Hb and deoxy-Hb signals with physiological artifacts) can be retained by using the inverse of the mixing matrix A^–1. Assuming that the physiological artifacts are statistically independent from the oxy-Hb and deoxy-Hb signals, they can be divided into separate ICs and removed from the fNIRS signal. For the decomposition of the fNIRS signal into ICs, the Second Order Blind Identification (SOBI) implementation of the ICA algorithm is used (Belouchrani et al., 1997).

Transfer function models

The use of TF models was introduced by Florian and Pfurtscheller (1997) in order to remove respiratory arrhythmia at movement induced changes of the HR. In NICA, the TF models are applied to remove the perturbations of the physiological artifacts from the oxy-Hb and deoxy-Hb signals. The models can be described with Eq 8:

$$X\left[n\right]=\sum _{u=0}^m{g}_{u}Y\left[n-u\right]+N\left[n\right](8)$$

X[n] describes the signal at each time point n, Y[n–u] the perturbation and N[n] the signal without influences. The parameters g_u can be estimated by minimizing the mean squared error E(N²[n]) (Florian and Pfurtscheller, 1997). After rearranging Eq 8, the signal without influences N[n] can be calculated with Eq 9:

$$N\left[n\right]=X\left[n\right]-\sum _{u=0}^m{g}_{u}Y\left[n-u\right](9)$$

Either one of these artifact correction methods can be chosen for the analysis (the choice may also depend on the measurement data and if physiological signals are available). However, according to Bauernfeind et al. (2014) the TF models have proven to be the most effective method.

Applying a Low-Pass Filter

In order to remove further artifacts (of technical or physiological origin), there is also the option to apply a low-pass filter (if the artifact lies in a frequency range above the frequency of interest). The used filter is a Butterworth filter of 8th order with a passband ripple of 3 dB and a stopband attenuation of 30 dB. This filter can be applied if there are any distortions visible in the frequency spectrum (e.g., originating from technical artifacts). The cut-off frequency is freely selectable, depending on the frequency range of the distortion. If physiological artifact correction is applied, the filter should not interfere with the physiological signals, which means. the cut-off frequency should not be in the frequency range of the physiological signals.

Grand Average and Region-of-Interest Analysis

With NICA it is possible to run a grand average analysis on the group level. Therefore, the oxy-Hb and deoxy-Hb signals of individual subjects are averaged over all subjects at each channel. Eq 10 shows the calculations for the oxy-Hb signals:

$$\mathrm{\Delta }{c}_{Hb}AVG=\frac{1}{N}\sum _{i=1}^N\mathrm{\Delta }{c}_{Hb}i(10)$$

The matrix with the average concentration changes of oxy-Hb Δc_HbAVG is calculated by summing up the individual concentration changes of oxy-Hb Δc_Hbi and then dividing them by the number of subjects N. The matrix notation is used, because the concentration change signals for all channels are calculated. The average concentration changes of deoxy-Hb are calculated in the same way.

At a grand average analysis, the oxy-Hb and deoxy-Hb signals from individual channels can be combined to a region-of-interest (ROI). The calculation of an ROI of the concentration changes of oxy-Hb is shown at Eq 11:

$$\mathrm{\Delta }{c}_{Hb}ROI=\frac{1}{C}\sum _{c=1}^C\mathrm{\Delta }{c}_{Hb}c(11)$$

The ROI of the concentration changes of oxy-Hb Δc_HbROI, as described at Eq 11, is calculated by summing up the concentration changes of oxy-Hb from individual channels Δc_Hbc and then dividing them by the number of channels C. The ROI of the concentration changes of deoxy-Hb is calculated in the same way.

Implementation

Overview and Requirements

NICA is available as a free and open MATLAB toolbox¹. The toolbox includes the MATLAB files to install and start the GUI, the files to perform the analysis steps, as well as a directory for user settings and a user manual. The toolbox was developed and tested with the MATLAB releases 2015b, 2017a and 2019b. There are also additional software packages included, which are necessary for certain data handling and signal processing steps: BioSig v1.5 (Schlögl and Brunner, 2008) and EEGLAB v2008beta (Delorme and Makeig, 2004).

As mentioned in section “FNIRS Measurement Device and File Requirements,” an HDR-file and an XDF-file are required for NICA to perform an analysis. The XDF-file contains all the recorded signals of a measurement: fNIRS, ECG, HR, BP, RF as well as the paradigm markers. The paradigm markers must be set before the measurements and sent to the LabRecorder via LSL-stream. Without the markers it is not possible to average the fNIRS signals over trials (epochs) or distinguish between different conditions.

The toolbox can be installed by running the included installation file. After installation, the directory of the toolbox is added to the MATLAB paths, and NICA can simply be started by typing nica in the MATLAB command window.

Surface of the GUI

As seen in Figure 1, the surface of the toolbox is divided into different panels: Information about the measurement data (A); settings for the pre-processing of the measurement data, correction methods for physiological artifacts, excluding noisy channels or trials, and options for visualizing the results of the analysis (B); information about the current analysis status (C); a push button to start the analysis (D); information about the individual analysis and signal processing steps (E); and menu bar items allow the user to load the measurement data and save the analysis settings (F).

FIGURE 1

Figure 1. The surface of the GUI is divided into different panels: Information about the measurement data (A); settings for the pre-processing of the measurement data, correction methods for physiological artifacts, excluding noisy channels or trials, and options for visualizing the results of the analysis (B); information about the current analysis status (C); a push button to start the analysis (D); information about the individual analysis and signal processing steps (E); and menu bar items allow the user to load the measurement data and save the analysis settings (F).

The functions of the panels and menu items are described at Table 1. Detailed information about the handling of the GUI will be provided in section “Performing an Analysis” based on an example analysis.

TABLE 1

Table 1. Description of the functions of the panels and menu items of NICA.

Performing an Analysis

In this section, the steps for performing an analysis of fNIRS data with NICA will be exemplified. The dataset for this example was recorded with the NIRScout 1624 measurement device and its corresponding recording software NIRStar. The experiment was divided into six runs; hence six HDR-files and six XDF-files were generated. With NICA, a single-run analysis, an analysis with several run-files combined, as well as a grand average analysis on the group level can be performed.

The participants had to perform either a motor execution (ME) task by applying force on a pressure sensor with the right hand, or a motor imagery (MI) task by imagining applying force on the pressure sensor. The force for both ME and MI tasks was instructed to be 20% (condition 1) or 40% (condition 2) of the participant’s maximum grip strength and had to be applied for 12 s. One trial consisted of a 5 s reference period (pre-task, looking at a reference cross), 12 s task time and 8 s resting period (looking at a black screen). In this example, a dataset with an ME task and 40% grip strength will be presented.

Selecting Analysis Path and Loading Measurement Data

Before an analysis can be started, the directory where the output-files will be saved has to be selected. After selecting the menu item Select Analysis Path at the menu bar entry File (Figure 1F), a dialog box opens where the preferred file path can be chosen. It is recommended, to use the directory of the raw XDF-files.

The raw XDF-files can be loaded under the menu bar entry Measurement Data (Figure 1F) with clicking on Load and then choosing HDR File and XDF File(s). It is possible to load all XDF-files from a recording session (e.g., when the experiment is divided into different runs) and to perform the analysis for all files at once. When the measurement files are loaded, the file names and file paths are displayed on the left side of the GUI (Figure 1A).

Defining Conditions

Another important step is to define the different conditions (tasks) of your measurement data. The condition names and markers can be set under the menu bar entry Measurement Data by clicking on the option Define Conditions. There, the number of conditions must be entered and a name and a marker value for each condition must be defined (in this case 20% for condition 1 and 40% for condition 2). In order to average the measurement data over all trials correctly, the condition marker must fit the paradigm marker (i.e., same value), which are stored in the XDF-file.

Adjusting Analysis Settings

The analysis settings can be adjusted on the right side of the GUI. The settings are divided into different panels, according to their functionality (Figure 1B). Under the menu bar entry User Settings (Figure 1F) it is also possible to load a file with predefined analysis settings.

Basics

The basic settings include the signal imaging method (to display the oxy-Hb and deoxy-Hb signals averaged over trials or continuously), selecting one of the predefined conditions, setting the correct number of trials (of the selected condition) and defining a task name. The selected condition and the task name will be included in the name of the generated output files.

One very important setting at this panel is the Probe Set. With this setting, the user can decide how many channels are displayed in the oxy- and deoxy-Hb figures (concentration change signals and topoplots). The number of channels selected should match the number of channels used for recording the fNIRS signals. The number of channels used for recording depends on the optode setup, more precisely on the number of sources and detectors used, as well as their spatial distribution. Figure 2 shows the different probe sets which can be selected with the number and the positions of the channels.

FIGURE 2

Figure 2. Different probe sets, which can be selected with the corresponding number and positions of the channels.

Signal processing

In this panel, the signal analysis method for removing the physiological artifacts must be chosen (CAR, ICA, or TF models). The selection can be to remove all artifacts (MTH-waves, respiratory patterns, and pulsation of the heart) or just one or two of them (depending on which physiological signals were recorded during the measurement). If the option Uncorrected is chosen, no physiological artifact correction will be applied. Furthermore, different filters can be applied by clicking on the corresponding check boxes: Removing baseline drift (high-pass filter), a notch filter (removing 50 Hz power line interference) and a low-pass filter with its cut-off frequency (in Hz, to remove technical artifacts or instead of a physiological artifact correction method).

Physiological artifacts

In this section, the search windows and the correction band for each physiological artifact (in the frequency spectrum) are defined. Typically, the MTH waves are in a frequency range from 0.07 to 0.13 Hz, the RF from 0.2 to 0.4 Hz and the HR lies between 1 and 1.5 Hz (Bauernfeind et al., 2014). It is very important, that the search windows and the correction bands are chosen correctly, otherwise the physiological signals cannot be found in the frequency spectrum and therefore cannot be corrected. It is recommended to use the standard values as mentioned above, and to look at the frequency spectrum after the analysis. If a physiological artifact was not removed correctly, the search windows and correction bands must be adapted.

Timing

Here, the trials can be divided (as chosen from the condition marker) into a pre-task, a task and a post-task period. The different periods are separated in the analysis results figures (oxy-Hb and deoxy-Hb signals) by black vertical lines. These settings only apply for the figures of the ox-Hb and deoxy-Hb signals; the averaging and of the trials is done according to the condition markers.

Artifacts

If there the data is noisy of artifact afflicted (artifacts of technical origin), single trials or channels can be manually excluded from the measurement data for the current analysis. If the box for excluding a trial is checked, an input dialog box opens where the trial number must be inserted. The excluded channels must be selected at the panel Channel after clicking on the corresponding checkbox. Selecting the option Consider Optode Failure opens an input dialog box, where the channel number of the faulty optode has to be typed in, as well as the channels with which the faulty optode will be interpolated.

Channels

At this panel, figures of the oxy-Hb and deoxy-Hb signals for individual channels can be created by simply clicking on the corresponding check box. There are also checkboxes to select all channels at once (All) and to delete the selection of the channels again (None). If the option Exclude Channels is selected, the channels to exclude can be chosen by the same way. With the radio buttons at Channel Options, it can be switched between settings of the channels to display and the channels to exclude.

Figure options

To visualize the analysis results, different types of figures can be generated. It is possible to create a figure with the physiological signals (biosignals figures) BP, RF, HR, and ECG, generate frequency spectrum figures (before and after artifact correction) and concentration change figures of the oxy-Hb and deoxy-Hb signals (optional including the standard deviation of the oxy-Hb and deoxy-Hb signals). Furthermore, a 2-D topographical map (topoplot) can be created, representing the changes of oxy-Hb and deoxy-Hb over the scalp. For that plot, the number of figures can be chosen, in order to show the topological concentration changes in specific time-points. The limits of the concentration range (in mM × mm) can also be set manually.

Starting the Analysis

If the measurement data is loaded, the analysis path selected and the analysis settings adjusted correctly, the text Ready to start Analysis will be displayed at the panel Analysis Status (Figure 1C). The analysis can then be started by clicking on the push button Start Analysis (Figure 1D). During the analysis, information about the individual analysis and signal processing steps will be provided at the panel Output Text (Figure 1E).

Error Handling

If an error occurs during an analysis, for example if a certain physiological signal necessary for the physiological artifact correction is not available, then the user gets informed about the error (Figure 3). The source of the error is displayed at the panel Analysis Status and the font color changes to red. More detailed information is displayed at the panel Output Text and at the popup window Error Dialog. In this example, the user is informed that there is no blood pressure data available, which would be displayed in the biosignals figure and used for the physiological artifact correction. In that way, the user knows how to avoid the error or how to adapt the analysis. Furthermore, the user is also informed if analysis settings are incomplete, incorrect (for example, a string instead of a value), or they do not match the dataset.

FIGURE 3

Figure 3. Error handling: The source of the error is displayed at the panel Analysis Status; more detailed information is displayed at the panel Output Text and at the popup window Error Dialog.

Grand Average Analysis

After all single subject analyses of a study are completed, it is possible to perform a grand average analysis of all the study subjects (group-level analysis). To apply this setting, the option Start at the menu bar entry Grand Average (Figure 1F) must be selected. After clicking on Start, the push button Start Analysis changes its name to Start Grand Average. Again, an analysis path for the new output files has to be defined. It is recommended to choose the same path as for the single subject analysis files.

If results of individual channels are preferred to be combined, it is possible to define regions of interest (ROIs). In order to do that, the number of ROIs must be defined and then the corresponding channels assigned to each ROI. This can be done by clicking on Define ROIs at the menu bar entry Grand Average.

Results

Description of the Output Files

During an analysis of fNIRS data with NICA, several output files and figures are created, which are saved automatically to the generated file path after the analysis. The folder structure is as follows: Name of the study/Subject ID/Run number/Condition

Depending on the selected figure options, figures with the physiological signals (BP, RF, HR, and ECG), frequency spectrum figures (all and/or single channels), concentration changes of oxy-Hb and deoxy-Hb (all and/or single channels) and a 2-D topographical map (topoplot), representing the changes of oxy-Hb and deoxy-Hb over the scalp, can be generated by the toolbox. The values of the concentration changes of oxy-Hb and deoxy-Hb, for each channel at each time point, are saved in a Microsoft Excel file, where they can be used for further statistical analysis. Furthermore, a MATLAB-file with all workspace variables is generated and the signal processing details of the analysis are saved in a text-file. All output files have the subject ID, the task name and the condition as prefix; the figure files have additionally the artifact correction method as suffix. Table 2 lists all output files which can be generated with NICA (without pre- and suffix).

TABLE 2

Table 2. Generated output files with file format and a description of their content.

Output Figures

In this section, the output figures generated with NICA are described on the basis of the example analysis of section “Performing an Analysis”.

The physiological signals are illustrated in Figure 4. For a better view of the details and the shape of the signals, the figure is zoomed in at a time-window of 20 s. The inspection of the figures can help to validate the values of the physiological signals. Figure 4A shows the BP waves with systolic (red line, BPsys) and diastolic (blue line, BPdia) values (in mm Hg). Both the diastolic BP values (>100 mm Hg) and the systolic BP values (>140 mm Hg) indicate high blood pressure (Hypertension). The movement of the thorax due to respiration is shown in Figure 4B and shows periodic extension and contraction of the respiration belt. The HR values (Figure 4C) lie between 60 and 80 bpm, which indicates a normal resting HR for adults. The patterns of the ECG signal with positive R-peaks and T-waves can be seen in Figure 4D.

FIGURE 4

Figure 4. Physiological signals: BP waves with systolic (red line, BPsys) and diastolic (blue line, BPdia) values (A); movement of the thorax due to respiration (B); HR values (C); and peaks of the ECG signal (D).

Figure 5 displays the power (in dB) of the frequency spectrum of the oxy-Hb (red line) and deoxy-Hb (blue line) signals, before (raw, thick lines) and after (clean, thin lines) artifact correction, averaged over all channels. The corrected frequency bands are highlighted in green. The peak of the MTH-waves between 0.05 and 0.2 Hz was detected and removed, as well as the peaks of the RF at around 0.3 Hz and the influence of the pulsation of the heart between 1 and 1.5 Hz.

FIGURE 5

Figure 5. Power of the frequency spectrum of oxy-Hb (red line) and deoxy-Hb (blue line) signals, before (thin line, raw) and after (thick line, clean) artifact correction, averaged over all channels. The corrected frequency bands are highlighted in green.

The concentration changes (in mM^∗mm) of oxy-Hb (red line) and deoxy-Hb (blue line) at channel 30, averaged over all trials can be seen in Figure 6. Channel 30 is located on the left hemisphere at the primary motor cortex (Brodmann area 4). The vertical lines at t = 0 and t = 12 indicate the start and the end of the right-hand ME task. During the task, there is an increase of the concentration of oxy-Hb accompanied by slight increase of the concentration of deoxy-Hb. Toward the end of the task, both concentrations of oxy-Hb and deoxy-Hb start decreasing again. These concentration changes show a typical behavior that indicates cortical activation (Strangman et al., 2002; Obrig and Villringer, 2003).

FIGURE 6

Figure 6. Concentration change of oxy-Hb (red) and deoxy-Hb (blue) at channel 30 (motor cortex, left hemisphere), averaged over all trials. The vertical lines at t = 0 and t = 12 indicate the start and the end of the right-hand ME task.

The 2-D topographical maps (topoplots), representing the changes of oxy-Hb (top row) and deoxy-Hb (bottom row) over the scalp, are shown in Figure 7. The concentration changes are divided into four time-windows (t1 = -5 to 0 s, t2 = 0–7 s, t3 = 7–14 s, t4 = 14–20 s). The colors orange/red indicate an increase of the concentration, green/blue a decrease. There is an increase of the concentration of oxy-Hb in the frontal and the parietal, motor related areas during the task, and only small changes in the concentration of deoxy-Hb. The highest increase in the concentration of oxy-Hb can be found at window t3, 7–14 s after task onset. These concentration changes corroborate the findings in Figure 6.

FIGURE 7

Figure 7. 2-D topographical maps, representing the changes of oxy-Hb (top row) and deoxy-Hb (bottom row) over the scalp, divided into four time-ranges (t1 = −5 to 0 s, t2 = 0–7 s, t3 = 7–14 s, t4 = 14–20 s). The colors orange/red indicate an increase of the concentration, green/blue a decrease.

In addition to the figures shown in this section, the frequency spectra and the concentration changes of oxy-Hb and deoxy-Hb of each individual channel can be generated as well.

Discussion

With NICA, a novel MATLAB toolbox for near-infrared spectroscopy calculations and analyses, it is possible to process and visualize fNIRS data in a straightforward and user-friendly way, as well as to apply physiological artifact correction methods. All output files and figures are saved automatically after an analysis in the generated analysis directory. Additionally, a proper error handling is implemented supporting an easy use of the toolbox. The surface of the GUI is divided into panels, which make it easy for the user to set and change different analysis settings that can also be saved and loaded for individual studies. Furthermore, NICA provides a comprehensive documentation of the analysis, which allows the user to follow and understand the individual signal processing steps.

For the correction of the physiological artifacts, three methods are implemented and described: CAR, ICA and TF models. According to Bauernfeind et al. (2014) using TF models with diastolic blood pressure as MTH-waves source is the most effective artifact correction method. However, it is still possible that in some cases the other methods can be more successful and lead to better results. In the example analysis in section “Performing an Analysis,” the best artifact correction method was TF models using the HR as MTH-waves source. It can be seen in Figure 5, that the physiological artifacts were successfully corrected. Nevertheless, it has to be mentioned that the artifact correction strongly depends on the signal quality of the physiological data. The correction of an artifact can only be as good as the recorded signal of this artifact.

The recorded data (fNIRS and physiological signals and paradigm markers) must be stored in an XDF-file, which is generated with the LSL recording software LabRecorder. It is possible to perform an analysis with fNIRS data exclusively (without additional physiological signals), but an artifact correction is only possible with the CAR method or with applying a low-pass filter then. Another important aspect is that the paradigm markers are set and recorded correctly. Without the markers it is not possible to divide the fNIRS signal into epochs and average it over trials or to distinguish between different conditions.

A limitation of NICA is the compatibility with different measurement systems. The software was developed for fNIRS data recorded with the NIRScout 1624 device. Although options for analyzing fNIRS data recorded with different measurement systems are available in NICA, they have not been tested or evaluated yet. Since the source code of the toolbox is freely available, it is possible to implement further or adapt existing functions.

Data Availability Statement

The datasets generated for this study can be found in the authors GitHub repository https://github.com/praggam/NICA.

Ethics Statement

Secondary data was used in this study from a preliminary study of Wriessnegger et al. (2017). This original study was approved by the local ethics committee (Medical University of Graz) and is in accordance with the ethical standards of the Declaration of Helsinki. After a detailed written and oral instruction, they gave informed written consent to participate in the study.

Author Contributions

PR and SW developed the structure and the content of the manuscript. PR wrote the text of the manuscript and made the figures. SW supervised the toolbox. GB wrote parts of the functions used in the toolbox. SW and GB proofread and corrected the manuscript. All authors read and approved the final manuscript.

Conflict of Interest

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

NICA was developed at the Institute of Neural Engineering at Graz University of Technology. It is a free and open-source software package that can be further developed and distributed non-commercially. The connection of the individual measurement signals (fNIRS, HR, BP, RF and paradigm markers) via LSL was installed by Dominik Bachmaier.

Footnotes

1. ^ https://github.com/praggam/NICA

References

Bauernfeind, G., Leeb, R., Wriessnegger, S. C., and Pfurtscheller, G. (2008). Development, set-up and first results for a one-channel near-infrared spectroscopy system / entwicklung, aufbau und vorläufige ergebnisse eines einkanal- nahinfrarot-spektroskopie-systems. Biomed. Tech. (Berl) 53, 36–43. doi: 10.1515/bmt.2008.005

CrossRef Full Text | Google Scholar

Bauernfeind, G., Wriessnegger, S. C., Daly, I., and Müller-Putz, G. R. (2014). Separating heart and brain: on the reduction of physiological noise from multichannel functional near-infrared spectroscopy (fNIRS) signals. J. Neural Eng. 11:056010. doi: 10.1088/1741-2560/11/5/056010

PubMed Abstract | CrossRef Full Text | Google Scholar

Belouchrani, A., Abed-Meraim, K., Cardoso, J. F., and Moulines, E. (1997). A blind source separation technique using second-order statistics. IEEE Trans. Signal Process. 45, 434–444. doi: 10.1109/78.554307

CrossRef Full Text | Google Scholar

Boas, D. A., Dale, A. M., and Franceschini, M. A. (2004). Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy. Neuroimage 23(Suppl. 1), S275–S288. doi: 10.1016/j.neuroimage.2004.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Coyle, S. M., Ward, T. E., and Markham, C. M. (2004). Physiological noise in near-infrared spectroscopy: implications for optical brain computer interfacing. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2004, 4540–4543. doi: 10.1109/IEMBS.2004.1404260

PubMed Abstract | CrossRef Full Text | Google Scholar

Coyle, S. M., Ward, T. E., and Markham, C. M. (2007). Brain-computer interface using a simplified functional near-infrared spectroscopy system. J. Neural Eng. 4, 219–226. doi: 10.1088/1741-2560/4/3/007

PubMed Abstract | CrossRef Full Text | Google Scholar

De Boer, R. W., Karemaker, J. M., and Strackee, J. (1986). On the spectral analysis of blood pressure variability. Am. J. Physiol. 251(Pt 2), H685–H687. doi: 10.1152/ajpheart.1986.251.3.h685

CrossRef Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Delpy, D. T., and Cope, M. (1997). Quantification in tissue near–infrared spectroscopy. Philos. Trans. R. Soc. Lond. B Biol. Sci. 352, 649–659. doi: 10.1098/rstb.1997.0046

CrossRef Full Text | Google Scholar

Duncan, A., Meek, J. H., Clemence, M., Elwell, C. E., Tyszczuk, L., Cope, M., et al. (1995). Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy. Phys. Med. Biol. 40, 295–304. doi: 10.1088/0031-9155/40/2/007

PubMed Abstract | CrossRef Full Text | Google Scholar

Elwell, C. E., Springett, R., Hillman, E., and Delpy, D. T. (1999). Oscillations in cerebral haemodynamics. Implications for functional activation studies. Adv. Exp. Med. Biol. 471, 57–65. doi: 10.1007/978-1-4615-4717-4_8

CrossRef Full Text | Google Scholar

Ferrari, M., and Quaresima, V. (2012). A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application. Neuroimage 63, 921–935. doi: 10.1016/j.neuroimage.2012.03.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Florian, G., and Pfurtscheller, G. (1997). Elimination von Atmungseffekten Auf Bewegungsinduzierte Änderungen Der Herzrate – elimination of respiratory effects on movement-induced cardiac response. Biomed. Tech. 42, 203–206. doi: 10.1515/bmte.1997.42.7-8.203

CrossRef Full Text | Google Scholar

Gratton, G., and Corballis, P. M. (1995). Removing the heart from the brain: compensation for the pulse artifact in the photon migration signal. Psychophysiology 32, 292–299. doi: 10.1111/j.1469-8986.1995.tb02958.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Huppert, T. J., Diamond, S. G., Franceschini, M. A., and Boas, D. A. (2009). HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain. Appl. Optics. 48, D280–D298. doi: 10.1364/ao.48.00d280

PubMed Abstract | CrossRef Full Text | Google Scholar

Jöbsis, F. F. (1977). Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters. Science 198, 1264–1267. doi: 10.1126/science.929199

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaur, M., and Singh, B. (2011). “Comparison of different approaches for removal of baseline wander from ECG signal,” in Proceedings of the International Conference & Workshop on Emerging Trends in Technology - ICWET ’11, ed. B. K. Mishra (New York, NY: Association for Computing Machinery), doi: 10.1145/1980022.1980307

CrossRef Full Text | Google Scholar

Koepchen, H. P. (1991). “Physiology of rhythms and control systems: an integrative approach,” in Rhythms in Physiological Systems. Springer Series in Synergetics, Vol. 55, eds H. Haken and H. P. Koepchen (Berlin: Springer), doi: 10.1007/978-3-642-76877-4_1

CrossRef Full Text | Google Scholar

Kothe, C. (2014). Lab Streaming Layer (LSL). GitHub. Available online at: https://github.com/sccn/labstreaminglayer (accessed February 17, 2020).

Google Scholar

Leff, D. R., Orihuela-Espina, F., Elwell, C. E., Athanasiou, T., Delpy, D. T., Darzi, A. W., et al. (2011). Assessment of the cerebral cortex during motor task behaviours in adults: a systematic review of functional near infrared spectroscopy (fNIRS) studies. Neuroimage 54, 2922–2936. doi: 10.1016/j.neuroimage.2010.10.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Leithner, C., and Royl, G. (2014). The oxygen paradox of neurovascular coupling. J. Cereb. Blood Flow Metab. 34, 19–29. doi: 10.1038/jcbfm.2013.181

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludwig, K. A., Miriani, R. M., Langhals, N. B., Joseph, M. D., Anderson, D. J., and Kipke, D. R. (2009). Using a common average reference to improve cortical neuron recordings from microelectrode arrays. J. Neurophysiol. 101, 1679–1689. doi: 10.1152/jn.90989.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Malonek, D., and Grinvald, A. (1996). Interactions between electrical activity and cortical microcirculation revealed by imaging spectroscopy: implications for functional brain mapping. Science 272, 551–554. doi: 10.1126/science.272.5261.551

PubMed Abstract | CrossRef Full Text | Google Scholar

Matthews, F., Pearlmutter, B. A., Wards, T. E., Soraghan, C., and Markham, C. (2008). Hemodynamics for brain-computer interfaces. IEEE Signal Process. Mag. 25, 87–94. doi: 10.1109/msp.2008.4408445

PubMed Abstract | CrossRef Full Text | Google Scholar

Obrig, H., Neufang, M., Wenzel, R., Kohl, M., Steinbrink, J., Einhäupl, K., et al. (2000). Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage 12, 623–639. doi: 10.1006/nimg.2000.0657

PubMed Abstract | CrossRef Full Text | Google Scholar

Obrig, H., and Villringer, A. (2003). Beyond the visible–imaging the human brain with light. J. Cereb. Blood Flow Metab. 23, 1–18. doi: 10.1097/01.WCB.0000043472.45775.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Okada, E., Firbank, M., Schweiger, M., Arridge, S. R., Cope, M., and Delpy, D. T. (1997). Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head. Appl. Optics 36, 21–31. doi: 10.1364/ao.36.000021

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Bauernfeind, G., Wriessnegger, S. C., and Neuper, C. (2010a). Focal frontal (de)oxyhemoglobin responses during simple arithmetic. Int. J. Psychophysiol. 76, 186–192. doi: 10.1016/j.ijpsycho.2010.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Ortner, R., Bauernfeind, G., Linortner, P., and Neuper, C. (2010b). Does conscious intention to perform a motor act depend on slow cardiovascular rhythms? Neurosci. Lett. 468, 46–50. doi: 10.1016/j.neulet.2009.10.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Power, S. D., Kushki, A., and Chau, T. (2011). Towards a system-paced near-infrared spectroscopy brain-computer interface: differentiating prefrontal activity due to mental arithmetic and mental singing from the no-control state. J. Neural Eng. 8:066004. doi: 10.1088/1741-2560/8/6/066004

PubMed Abstract | CrossRef Full Text | Google Scholar

Santosa, H., Zhai, X., Fishburn, F., and Huppert, T. (2018). The NIRS brain AnalyzIR toolbox. Algorithms 11:73. doi: 10.3390/a11050073

CrossRef Full Text | Google Scholar

Schlögl, A., and Brunner, C. (2008). BioSig: a free and open source software library for BCI research. Computer 41, 44–50. doi: 10.1109/mc.2008.407

CrossRef Full Text | Google Scholar

Scholkmann, F., Kleiser, S., Metz, A. J., Zimmermann, R., Pavia, J. M., Wolf, U., et al. (2014). A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology. Neuroimage 85(Pt 1), 6–27. doi: 10.1016/j.neuroimage.2013.05.004

CrossRef Full Text | Google Scholar

Soraghan, C., Matthews, F., Markham, C., Pearlmutter, B. A., O’Neill, R., and Ward, T. E. (2008). A 12-channel, real-time near-infrared spectroscopy instrument for brain-computer interface applications. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2008, 5648–5651. doi: 10.1109/iembs.2008.4650495

PubMed Abstract | CrossRef Full Text | Google Scholar

Strangman, G., Culver, J. P., Thompson, J. H., and Boas, D. A. (2002). A quantitative comparison of simultaneous BOLD fMRI and NIRS recordings during functional brain activation. Neuroimage 17, 719–731. doi: 10.1006/nimg.2002.1227

PubMed Abstract | CrossRef Full Text | Google Scholar

Tachtsidis, I., and Scholkmann, F. (2016). False positives and false negatives in functional near-infrared spectroscopy: issues, challenges, and the way forward. Neurophotonics 3:031405. doi: 10.1117/1.nph.3.3.030401

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, T., Takikawa, Y., Kawagoe, R., Shibuya, S., Iwano, T., and Kitazawa, S. (2011). Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task. Neuroimage 57, 991–1002. doi: 10.1016/j.neuroimage.2011.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Villringer, A., Planck, J., Hock, C., Schleinkofer, L., and Dirnagl, U. (1993). Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults. Neurosci. Lett. 154, 101–104. doi: 10.1016/0304-3940(93)90181-j

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, T., Fujiwara, T., and Suzuki, S. (2018). “NIRS-SPM analysis of body schema modification and performance of body motion,” in Proceedings of the 2018 11th International Conference on Human System Interaction (HSI) (Gdansk: IEEE), doi: 10.1109/hsi.2018.8430778

CrossRef Full Text | Google Scholar

Wolf, M., Ferrari, M., and Quaresima, V. (2007). Progress of near-infrared spectroscopy and topography for brain and muscle clinical applications. J. Biomed. Optics 12:062104. doi: 10.1117/1.2804899

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, M., Wolf, U., Toronov, V., Michalos, A., Paunescu, L. A., Choi, J. H., et al. (2002). Different time evolution of oxyhemoglobin and deoxyhemoglobin concentration changes in the visual and motor cortices during functional stimulation: a near-infrared spectroscopy study. Neuroimage 16(Pt 1), 704–712. doi: 10.1006/nimg.2002.1128

PubMed Abstract | CrossRef Full Text | Google Scholar

Wriessnegger, S. C., Kirchmeyr, D., Bauernfeind, G., and Müller-Putz, G. R. (2017). Force related hemodynamic responses during execution and imagery of a hand grip task: a functional near infrared spectroscopy study. Brain Cogn. 117, 108–116. doi: 10.1016/j.bandc.2017.06.010

CrossRef Full Text | Google Scholar

Zhang, Q., Brown, E. N., and Strangman, G. E. (2007). Adaptive filtering to reduce global interference in evoked brain activity detection: a human subject case study. J. Biomed. Optics 12:064009. doi: 10.1117/1.2804706

PubMed Abstract | CrossRef Full Text | Google Scholar