Tutorial on platform for optical topography analysis tools

[+] Author Affiliations
Stephanie Sutoko, Hiroki Sato, Atsushi Maki, Masashi Kiguchi, Yukiko Hirabayashi, Hirokazu Atsumori, Akiko Obata, Tsukasa Funane, Takusige Katura

Hitachi Ltd., Research and Development Group, 2520 Akanuma, Hatoyama, Saitama 350-0395, Japan

Neurophoton. 3(1), 010801 (Jan 11, 2016). doi:10.1117/1.NPh.3.1.010801
History: Received June 16, 2015; Accepted December 2, 2015
Text Size: A A A

Open Access Open Access

Abstract.  Optical topography/functional near-infrared spectroscopy (OT/fNIRS) is a functional imaging technique that noninvasively measures cerebral hemoglobin concentration changes caused by neural activities. The fNIRS method has been extensively implemented to understand the brain activity in many applications, such as neurodisorder diagnosis and treatment, cognitive psychology, and psychiatric status evaluation. To assist users in analyzing fNIRS data with various application purposes, we developed a software called platform for optical topography analysis tools (POTATo). We explain how to handle and analyze fNIRS data in the POTATo package and systematically describe domain preparation, temporal preprocessing, functional signal extraction, statistical analysis, and data/result visualization for a practical example of working memory tasks. This example is expected to give clear insight in analyzing data using POTATo. The results specifically show the activated dorsolateral prefrontal cortex is consistent with previous studies. This emphasizes analysis robustness, which is required for validating decent preprocessing and functional signal interpretation. POTATo also provides a self-developed plug-in feature allowing users to create their own functions and incorporate them with established POTATo functions. With this feature, we continuously encourage users to improve fNIRS analysis methods. We also address the complications and resolving opportunities in signal analysis.

Optical topography/functional near-infrared spectroscopy (OT/fNIRS) noninvasively measures cerebral hemoglobin concentration changes using two or more wavelengths (650 to 950 nm), where oxygenated- and deoxygenated-hemoglobin (O2Hb and HHb) have different absorption coefficients.14 Related to cerebral energy metabolism,5 oxygen demand and consumption increase while the brain is activated. Therefore, O2Hb and HHb concentration changes (ΔCO2Hb and ΔCHHb) may indicate brain activation. The fNIRS has been used for more than 20 years and potentially has broad applications such as for neurologic and psychiatric disorders, cognitive/developmental psychology, and rehabilitation68 due to advances in safety, convenience, cost-effectiveness, and temporal resolution compared with functional magnetic resonance imaging (fMRI).4,912 In spite of these advantages, fNIRS has limited spatial resolution due to scattered light in the superficial brain layers (scalp and skull), and the signals are susceptible to motion and systemic noise.1316 Regarding the noise problem, the main issue lies with bias interpretation of functional signals. Though many methods, such as evoked-signal modeling,17,18 coherent systemic separation,1921 phase synchronization,2224 and hemodynamic modalities,25,26 have been proposed to extract the brain activation signals using various available platforms and analysis tools,18,2731 the suitability of an analysis method that may depend on the performed tasks and participant has not been extensively examined. For example, signal modeling might be different depending on the task or systemic, and the noise separation method will affect the analysis of resting-state and motoric-task (e.g., exercise) studies differently.

We developed a MATLAB®-based software called platform for optical topography analysis tools (POTATo), an integrated data analysis environment that assists users in comprehensively analyzing fNIRS data for multiple purposes.3239 POTATo was established by integrating four essential element functions, the compatibility of data handling, ease of exploratory data analysis (EDA) and model diagnostics (MD), and result/data presentation in a holistic platform, as illustrated in Fig. 1. These functions enable users to (1) organize particular datasets in separate POTATo projects, (2) modify the projects by importing additional data, exporting, or merging with other projects, (3) conduct EDA to roughly understand raw data,40 (4) visualize raw data, transient-preprocessed signal, and statistical results, and (5) verify MD via statistical analysis functions presenting parallelism between EDA and MD. Engaged by many researchers around the world, fNIRS analysis is being continually improved and will result in novel and prospective analysis methods. Hence, POTATo is provided with a self-developed plug-in feature allowing users to create functions, install them in the POTATo package, and either use them together or compare them with other POTATo functions.

Graphic Jump Location
Fig. 1
F1 :

Concept of platform for optical topography analysis tools (POTATo) software with integrated analysis functions.

Compared with widely used analysis software (e.g., HomER2 and NIRS-SPM),27,28 POTATo does not have a definite procedure for analyzing data, which provides (1) higher flexibility in selecting and arranging preprocessing functions, (2) accessibility for users to directly involve in creating functions and evaluating processed data/results, (3) an integrated analysis platform to simultaneously preprocess, analyze, and infer a dataset, (4) analysis transparency enabling users to export processed data/results, and (5) feasibility of incorporated analysis with other domains (e.g., statistics software). Therefore, we anticipate that POTATo can convey insightful perspectives and strategies in data analysis and effectively promote further fNIRS research. In this tutorial, we describe the conceptual architecture of POTATo, guidelines in operating and exploring the POTATo graphical user interface (GUI) and functions, and an example of fNIRS data analysis using an actual cognitive study of working memory (WM) tasks.

POTATo is an analysis platform that invites users to conduct data analysis by either selecting or creating preprocessing and analysis methods. However, this requires users to have experience in handling and analyzing fNIRS data; therefore, POTATo promotes two modes, Normal and Research, having different features, flexibility, and purposes. In Normal mode, functions are compiled by orderly arranged preprocessing analysis and instantly applied to a dataset without any further adjustment. In contrast, Research mode enables users to design sequential preprocessing and analysis functions. This mode also allows users to specify the parameters and easily observe processed data transformation while determining the functions. Research mode is more useful for researchers who would like to delve into and verify hypotheses, whereas Normal mode is more suitable for automatic analysis of studies that have been previously optimized such as depression level diagnosis.32,34,37,39 Therefore, the optimum analysis function can be directly implemented to a dataset with the same experimental paradigm and objective. This tutorial presents an understandable explanation of fNIRS data analysis including how to preliminarily assess the raw data, choose the necessary preprocessing functions, adjust parameters, and infer the conclusion; thus we encourage users to practice and analyze the data in Research mode.

Research mode organizes data analysis in three structural interfaces starting from signal preprocessing, data grouping, and statistical testing. Figure 2 illustrates the Research mode GUI together with data flow. The Preprocess interface permits users to arrange a list of preprocessing functions called a recipe for each single datum. Moreover, the layout of raw and processed data can be evaluated through visual examination. The transient data during preprocessing leads to a fundamental basis for evaluating the effectiveness of signal preprocessing in minimizing noise. In this step, users will get used to observing and apprehending information from raw data (i.e., EDA) and processed data, and specify the recipe. Since completing the preprocessing step, the following step is done in the Summary Statistic Computation interface to compartmentalize data, determine the analysis method (e.g., averaged cerebral hemoglobin concentration changes within predetermined interval), analyze the processed data according to the selected method, and form the grouped data summary for preparing the statistical analysis. The grouped data summary will be analyzed further in the Statistical Test interface, with which users are able to carry out the analysis by conveniently controlling the statistic parameters. The statistical analysis results are supposed to guide users to a conclusion of proving or rejecting the primary hypotheses, i.e., MD.

Graphic Jump Location
Fig. 2
F2 :

Sequential steps in functional near-infrared spectroscopy (fNIRS) data analysis using POTATo software including preprocessing, creating group analysis, and statistical analysis in three different interfaces.

Figure 3 shows how POTATo works on raw data to further understand data processing. After importing a dataset, POTATo stores the original Raw data (MAT-files) and constructs an empty Recipe (MAT-files) for each datum. Instead of modifying the original Raw data, a user solely sets, edits, and overwrites the Recipe. POTATo loads the Raw data and processes according to the function list in the Recipe. Any deletion done in POTATo (menu Edit → Data Delete) merely removes recipes. It, therefore, will not eliminate the Raw data from the POTATo project. Considering how important the Raw data is, this structural arrangement is selected to avoid the mistake of permanently dismissing the Raw data. If a user unintentionally removes the Recipe, POTATo can track down the missing Recipe by comparing the list of imported Raw data with the currently remaining Recipe. With the user’s permission, POTATo will create a new Recipe. However, the modification on probe location and stimulus onset time information may be permanently executed on the Raw data. The user needs to be careful in performing this modification. In detail, the description of POTATo features, functions, and stepwise analysis procedure are covered in Sec. 3.

Graphic Jump Location
Fig. 3
F3 :

POTATo processing in modifiable Recipe files and safely stored Raw data.

POTATo can be downloaded (available in Ref. 41) for free. This tutorial is on POTATo version 3.8.0 beta, and we are continuously improving POTATo features. The POTATo package is wrapped in a .zip file called P3_files_EN. This includes two POTATo versions depending on the user’s MATLAB® version and a brief user’s guide for installation and POTATo processing. To get started, the POTATo package (P38 folder) is set on the MATLAB® path and a user can launch the platform by calling function ‘P3’ in the command window. For first installment, the platform default will show the selectable window to either load an existing project or create a new project. Beginning with creating a new project, the user is asked to name the project and import a dataset collected from an fNIRS instrument (e.g., ETG and WOT series, Hitachi Medical Corporation and Hitachi High-Technologies in .csv file; see Appendix 1 for importable data format details) or exported from other POTATo projects (Project → Modify → Export in .mat file). Once imported, measurement information, such as channel number, position configuration, and sampling rate (recognizable for ETG and WOT series by default) is transferred and the channel-wise intensity signals (e.g., analog-to-digital converted signals from optical detected intensity for each wavelength) are converted following a modified Lambert Beer equation1,42 into POTATo raw data containing ΔCO2Hb, ΔCHHb, and ΔCHb-total data. Therefore, POTATo mainly processes cerebral hemoglobin data instead of intensity data.

After completing data import, the Research mode GUI will come up by default, as depicted in Fig. 4. For creating, accessing other projects, modifying the current project, or closing the POTATo platform, the user can use the main menu (Project → New project/Open/Modify/Close/Exit). To change the POTATo mode, edit the list of frequently used functions called My List, and set the project directory, the menu Setting can be used to execute these functions. On the left side of the GUI, there are several items such as Search text-field, Data List, Data Information, and Status columns, Extended button [Fig. 4(A)], and Position Setting, Mark Setting, and Multifunctional buttons [Fig. 4(B)]. The user can easily find data in the current project by typing the file name in the Search field. Moreover, by clicking the Extended button, another window called Extended Search [Fig. 4(C)] is launched. From that window, the user can choose data characteristics (e.g., gender, age, and measurement date) from the Key search list, set the parameter option, and click the Add button. More than one characteristic can be managed to define the search option together with the relationship among those characteristics; for example, female participants who are older than 20 years old (AND logical operation) or either female participants or participants who are older than 20 years old (OR logical operation). After selecting the characteristics, the Data List on the Preprocess interface will consequently change and the order of the Data List can be sorted in ascending or descending order. Modification of the Key search list can be done by exporting the current Key search file (.csv file), adjusting the listed characteristics, and reimporting the Key search file.

Graphic Jump Location
Fig. 4
F4 :

Preprocess interface of POTATo Research Mode consists of search feature with (A) Extended button, (B) Multifunctional button, (C) Extended Search GUI that appears by clicking Extended button, (D) analysis interface guide, (E) preprocessing function directory, (F) list of functions in particular directory, (G) Recipe editing button, (H) visualization tool, and (I) data exporting tool.

Returning to the main GUI of Research mode, the Data List, Data Information, and Status boxes contain the list of Recipes in a project, information of a particular datum, and status of POTATo operation, respectively. Even though a user is prevented from permanently deleting the POTATo Raw data from the POTATo GUI, basic information, such as channel position and onset stimulus time, can be irreversibly altered in Position and Mark Setting. Measurement systems are compiled together with probe information, and POTATo adapts to the channel position (half the distance between source and detector by default) in presenting the mapping results. However, it can be flexibly readjusted depending on the targeted brain region and participant. If studies on adults use 47 channels to thoroughly cover the prefrontal cortex (PFC), those on an infant may necessitate less than half that number. A user needs to modify the channel position from the default position. The Position Setting button enables a user to rearrange and customize the channel position. Similar to the Position Setting, Mark Setting enables a user to permanently remove and modify task-related stimulus time of the POTATo Raw data. Task-related analysis mainly focuses within the stimulus period; thus, the marker details are important in data analysis. In addition, marker modification is globally valid for all channels. Because of fixed alteration in the POTATo Raw data, a user is expected to be more careful in conducting these functions. Another button called Multifunctional presents several functions allowing the user to deliberately select the intended functions by pointing a mouse cursor on the button and right-clicking. The specific functions (e.g., spatially-registered position, project repair, and plug-in installment tool) are explained in Sec. 3.1.6.

Initial Exploratory Data Analysis and Temporal Preprocessing

As described in Sec. 2, Research mode includes three steps, Preprocess, Summary Statistic Computation, and Statistical Test, which are constructed in separate interfaces. There are three buttons [Fig. 4(D)] on the upper right of the GUI indicating each step to direct toward particular functional interfaces. To demonstrate fNIRS data analysis in POTATo, this tutorial presents practical spatial and verbal WM studies conducted by Sato et al.34 and Aoki et al.37 The fNIRS system (ETG-7100 by Hitachi Medical Corporation, Japan) was equipped for measurement with a 3×10 sensor-detector arrangement (47 channels). Detailed information of related tasks can be referred to in these studies and data from participants were obtained according to the regulations of the internal review board of the Central Research Laboratory, Hitachi, Ltd., following receipt of written informed consent. The total amount of analyzed data from 160 participants was preliminarily screened using the Edinburgh Handedness Inventory;43 five participants were not included because of left handedness and ambidexterity. To begin, EDA was conducted through blunt-visual observation in a singular datum and averaged data. Ignoring the noise probability and familiar hemodynamic response function (HRF), insight into raw data would at least allow the understanding of data tendency and give an idea or strategy to analyze those data.

The data are listed and classified according to the type of WM tasks by practicing the Extended Search feature. The Plot average button will appear once multiple data are selected. This provides the channel-wise plot and map of average brain activation from all selected subjects. To understand brain activity during performed tasks, the concerned time frame was limited to the stimulus period that also includes pre- and poststimulus of 1 and 16 s, respectively, in these WM tasks. During one experiment, there were eight stimuli provided for 8.5 s and the POTATo function called Blocking accommodated the extraction of the task-related interval and represented it as the average of cerebral hemoglobin concentration change over several trials. While still dragging all spatial and verbal WM data separately, the Blocking function is selected from the list. When the Add button is clicked, the parameter option appears and the predetermined periods of pre- and poststimulus are set. The Plot average button is once again clicked resulting in spatial-evoked signals in 25.5 s, as illustrated in Fig. 5(a). Regardless of the task, both waveform and temporal preprocessing effects were similarly presented (verbal-evoked task not shown).

Graphic Jump Location
Fig. 5
F5 :

(a) Relative concentration change (eight averaged trials) in oxygenated hemoglobin (oxy-Hb/O2Hb, red-lines) and deoxygenated hemoglobin (deoxy-Hb/HHb, blue-lines) in temporal trial span without any preprocessing steps and (b) temporal processed signal for spatial WM task from all subjects.

In the POTATo plot GUI, a user can control the range of x- and y-axes that describe time and hemoglobin concentration changes. On the upper left of the GUI, there is a list of available signal types including ΔCO2Hb, ΔCHHb, and ΔCHb-total. This allows the user to select which signals are necessarily presented in a plot. For clear visualization, the plots in this tutorial mainly consist of ΔCO2Hb (red line by default) and ΔCHHb (blue line by default). According to Fig. 5(a), several channels, mainly located in the upper PFC, exhibited heavy noise with high amplitude. Other than those channels, channels in the center to lower parts clearly delivered incremental ΔCO2Hb with a particular two-peaks pattern while executing both spatial and verbal WM tasks. This increase generally corresponds to brain activation; however, we cannot straightforwardly determine particular regions of interest (ROIs) and eliminate other potential channels because of intrasubject variability in systemic activity.4446 A dissimilar baseline among channels also occurred and could result in misleading information of brain activation. Therefore, comparing brain activity with map and plot visualizations is insufficient to infer any conclusions from spatial and verbal WM tasks. We investigated further by using the preliminarily obtained tendency from mean-brain activation information and conducting statistical significance analysis. The statistical results will help us conclude these hypotheses.

Before moving on to statistical issues, there are points that need to be addressed. Regarding noise occurrence, no information could be derived from the upper channel, possibly hiding the functional brain signal. In spite of noticeably high-amplitude noise, integrated noise with a smaller amplitude and short interval might affect the prone channels. For example, the systemic activation, overlaid cortical absorption, and task-related motion (e.g., tilting head up, enforced respiration, jaw motion) deliver random interferences that are hardly differentiated and separated from brain activity.4752 The fluctuating interferences are generally diverse among people and the negligence of this noise will convey high variability and inconsistency. Therefore, noise removal can essentially improve data reliability. With acceptable methods minimizing signal loss, noise elimination should be attempted. Further preprocessing, as explained later, was done; resulting in reduced-noise signals, as illustrated in Fig. 5(b).

As another point regarding strategic data analysis, we would like to introduce two approaches, average sample- and model-based analyses for proceeding with POTATo. In average sample-based analysis, the activation is interpreted from the significantly increased average amplitude of activation-related period compared with that of baseline period (i.e., prestimulus) without any predetermined assumption on a signal waveform. Despite its simplicity, the temporal information outside a predetermined activation-related period is neglected. On the contrary, statistical parametric mapping (SPM), a well-known spatiotemporal normalization mapping based on mass univariate analysis across voxels in positron emission tomography5355 and fMRI,17,56,57 is used to evaluate the correlation between signal and estimated HRF. The assumption of complex HRF convolution or the dynamic cerebral autoregulation concept (e.g., myogenic, metabolic, and neurogenic activities)5861 differentiates this SPM method with an average sample-based approach. The correlation representing the signal-model likelihood is assessed based on a linear model recognized as the general linear model (GLM).6264 This model accommodates the voxel-wise analysis by fitting the response variables, yj (j=1,2,,J observed time points) with explanatory variables, xij (i=1,2,,I regressor models), weighted parameters, βi corresponding to each regressor, and independent error, ϵj as yJ=xIJβI+ϵJ.6568 Unlike the conventional time-invariant model in early SPM, the temporal adaptive (i.e., delayed and dispersed evoked response) emphasizes the advantages of focal ROIs, task variability, subject and signal (e.g., ΔCO2Hb, ΔCHHb) characteristics in the GLM model.6971 By evaluating the model resemblance, the irrelevant signal (e.g., physiological low-frequency oscillation) and prompt high amplitude due to motion artifacts reduce statistical robustness.72 This also emphasizes the importance of accomplishing appropriate preprocessing.

To conduct signal preprocessing, we now discuss the Preprocess interface. The right window of the Preprocess interface provides a tool for preprocessing functions. Several plug-ins have already been installed in POTATo and listed in the Recipe editor. The upper scroll-box is the directory function and set as My List by default [Fig. 4(E)]. The function list of a related directory is then accessed through the lower scroll-box [Fig. 4(F)]. As described above, the function list on My List can be modified in the menu (Setting → My List). A user can also explore more functions by changing the function directory to All Filter. By simply choosing one of the functions on the lower scroll-box and clicking the Add button, POTATo will insert that function into the Recipe, as shown in the lower input-box. Additional functions will be listed and executed sequentially. In the application of signal filtering and fitting functions, the order becomes important and the sequence modification is readily provided by selecting the particular function on the Recipe and clicking the UP-DOWN buttons [Fig. 4(G)] on the left side of the scroll-boxes.

To optimize the preprocessing step, a user can use a trial-and-error method in designing the preprocessing recipe. The optimized recipe is determined not only by selecting the necessary functions but also by setting the proper function parameters. There are three functional buttons, Change, Enable/Disable, and Remove [Fig. 4(G)], for altering the function parameters, activating/deactivating the functions, and removing the functions from the Recipe, respectively. Attempting optimization requires visual observation for comparing the signal transformation and reconstruction before and after implementing the functions. The POTATo data visualization plug-in [Fig. 4(H)] is located on the lower-left of the GUI, which provides a scroll-box consisting of available plotting types, such as channel check, line plot, activation map, and statistical result map (h-value). These plots are not always applicable depending on the preprocessing functions. If the statistical test function is not included in the Recipe, the statistical result map will not work. The Edit button next to the plotting scroll-box allows a user to modify the plot function by editing, adding, or even removing the GUI elements. By clicking the Draw button, POTATo analyzes the signal according to the Recipe and presents the processed signal according to the selected plot type. For a user who wants to investigate further preprocessing functions and data, the M-file script and Export Workspace (WS) buttons [Fig. 4(I)] are provided for exporting the code into the MATLAB® script file and the processed data into the MATLAB® variables, respectively. Setting the visual plot and clicking the Draw check-box beside the M-file script button includes the plotting function in the Recipe code script. After becoming familiar with the Preprocessing interface, the following subsections will explain the basic selection of preprocessing functions and those applications to WM tasks.

Noise detection and rejection

Despite highly supervised and controlled experiments, motion noise, such as unconscious head tilting and frowning, might occur. To examine the occurrence of motion artifact in a signal, the frequency spectrum is assessed beforehand via the Spectrum Check layout [plotting scroll-box → Spectrum Check → Draw button; Fig. 4(H)]. The spectrum check lies on the discrete-time Fourier transform decomposing the temporal signals and congregating frequency-based components. The spiky spectrum toward the amplitude of ΔCO2Hb/HHb at an unusually high frequency might be correlated with motion noise. In addition, the motion artifact is presumably bounded with wider affected areas. If the respective spike also occurs at the same time in adjacent channels, it might coherently support the forejudged assumption of the motion artifact. Figure 6(a) illustrates the Spectrum Check GUI of the ΔCO2Hb signal (channel 11, 3×10 S-D arrangement) from a representative subject performing a spatial WM task. Similar to other plotting GUIs, several elements, such as channel, available signal list-box, and axis controller, are provided [Figs. 6(A)6(C)]. According to the ΔCO2Hb plot (lower plot), which showed a sharp peak in the third poststimulus interval (94 s), there was also a spike (shown with an arrow) in a frequency over 2 Hz observed in ΔCO2Hb amplitude spectrum (upper plot) that might be related to the ΔCO2Hb spike. We can closely examine the spike by adjusting the time-frame (88 to 104 s) using the time scroll-bar at the bottom the GUI [Fig. 6(D)], as shown in Fig. 6(b). The neighboring channels (e.g., channels 21 and 31) also showed similar phenomena at that time point [Figs. 6(c)6(e)]. Therefore, that spike was likely due to a motion artifact.

Graphic Jump Location
Fig. 6
F6 :

(a) Spectrum Check layout in representative ΔCO2Hb signal from channel 11 while performing spatial WM task. Upper plot shows ΔCO2Hb amplitude spectrum versus frequency and lower plot illustrates ΔCO2Hb versus time. This allows user to (A) set displayed channel, (B) one or more signal types, (C) x- and y-axes, and (D) displayed time interval. Suspected spiky noise was observed in poststimulus of 3rd task trial (shown with an arrow). (b) Detail observation of high-frequency noise around 94 s at channel 11 and its adjacent channels (c) 21 and (d) 31 as (e) diagonally located.

While visual evaluation is important to thoroughly observe signals, it is merely practical for small-scale data analysis; otherwise, it will be time consuming and very laborious. Following Pena’s method using signal threshold and spike recognition,73 a POTATo function called Motion Check can identify and mark the determined noise depending on the adjusted parameters. Once added into the Recipe, the function GUI (Fig. 7) appears and the user is asked to specify several parameters including filtering method [e.g., ButterWorth, bandpass2/Fast Fourier Transform (FFT), or none], bandpass frequency of filtering, signal threshold, concerned data-kind, and interval of spike recognition. The term “data-kind” is frequently used in other function GUIs and structural data, which numbers 1, 2, and 3 denote ΔCO2Hb, ΔCHHb, and ΔCHb-total, respectively.

Graphic Jump Location
Fig. 7
F7 :

Motion Check interface conveniently confirms noise occurring at particular time points for each channel and data-kind (1 for ΔCO2Hb, 2 for ΔCHHb, and 3 for ΔCHb-total). Impact on channel can be immediately monitored as parameter justification. Detected ΔCO2Hb noise is denoted with cross-marks on channels 11 and 21.

As discussed above, the bias interpretation of the brain activation signals might be a limitation of fNIRS data analysis. This can be caused by unrelated activation signals convoluted together. In a simple method, the signal convolution problem is usually resolved by applying the filter to minimize the effect of particular frequency-based components toward the entire time-course signal. The band-passes are set to 0.02 and 1 Hz to reduce noise related to low-frequency oscillation and heartbeat effects, respectively. If the user inputs 0.1  mM·mm, kind 1, and 2 points in the threshold, data-kind, and interval point boxes, respectively, POTATo will look into the ButterWorth-filtered ΔCO2Hb signal for each channel, detect whether there are any ΔCO2Hb deviations over the 0.1-mM·mm threshold within two consecutive sampling points, and mark those exceeded points as noise. These parameters can be customized depending on the data condition. The marked flaw signal within a stimulus-interest interval will be excluded in further analysis. Either extremely strict or lenient motion detection will result in no signal being eligible for analysis or no noise being rejected, and both phenomena should be avoided. The user may select “None” for evaluating the POTATo Raw data without any filter; however, this can overestimate noise and ineffectively eliminate data.

On the right side of the GUI, two plots are provided for displaying the POTATo Raw data (bold-colored line), filtered signal (thin-colored line), and marked noise (cross-mark) from two selectable channels. As illustrated in Fig. 7, the line plots of channels 11 and 21 were presented. In addition, the spiky noise at 94 s, several time-points marked with cross-marks were mutually detected at an over-intensified amplitude (over 0.1  mM·mm). These results also improved this function’s efficiency compared with subjective visual observation. After finishing the parameter setting, the OK button is clicked to save the set parameters and this function is added to the Recipe. Since applying this function does not afford signal reconstruction and the filtering method solely works for initial signal smoothing before noise detection, another function called Band Filter is required.

Signal Filtering

As depicted in Fig. 8, the Band Filter function has several parameter inputs such as filter type (e.g., FFT and ButterWorth), band-filtering (e.g., high-pass, low-pass, and both filters), and filter dimension,18 especially for the ButterWorth filter on the left side of the GUI. The frequency band criteria are then specified in the input-boxes of the high-pass filter and low-pass filter. Two scroll-boxes below the low-pass filter input-box allow the user to select and visualize the channel-wise signals (e.g., ΔCO2Hb, ΔCHHb, and ΔCHb-total) that are outlined on the right side of the GUI. There are four plots including power density [upper two plots, Figs. 8(A) and 8(B)] and ΔCHb plots [lower two plots, Figs. 8(C) and 8(D)] for raw [Figs. 8(A) and 8(C)] and filtered signals [Figs. 8(B) and 8(D)] in a selected channel. These plots are expected to facilitate direct monitoring while justifying and customizing parameters.

Graphic Jump Location
Fig. 8
F8 :

Band Pass function GUI includes several parameters such as filter type, frequency, displayed channel, and signal types. Signals of channel 21 were observed in (A) and (B) spectrum plots and time course ΔCHb (C) and (D) reconstruction plots for Raw data (A) and (C) and filtered signals (B) and (D).

Task-related signal

POTATo was designed by adopting the block design analysis74,75 and several plug-ins particularly process continuous, block, or both data. Continuous data is considered as time-course data during an experiment, whereas block data is extracted data from continuous data in accordance with the stimulus markers. Because there are two types of markers, block and event markers, a user should confirm the task design beforehand. A block marker denotes the continuous task within a specific interval, as done in WM studies, while an event marker describes the performed task as a single action such as one-press button. The time period for single block data called an epoch is predetermined and composed of pre-, task-, and poststimulus intervals. For these WM tasks in which the stimulus period lasted for 8.5 s, pre- and poststimulus intervals were set to 1 and 16 s, respectively. The Blocking function is added to create block data. In addition, the interval inputs (e.g., pre- and poststimulation), a user also needs to specify the marker number. When the task paradigm is designed consisting of more than one type of task, POTATo recognizes the different stimulus types by distinguishing the marker number. Regardless of the marker number, the user can input All for this parameter.

To visualize the epoch, we can select either Channel Check or Line Plot layouts. Because the Motion Check function has already been applied, the displayed signal comes from the average of free-noise epochs. To observe the number of free-noise epochs, the Line plot (block overlay) layout illustrates the channel-wise visualization with the number of line-plots, which each line-plot defines a single free-noise epoch, as depicted in Fig. 9. The empty channel plot indicates that all epochs are affected by noise. If channels mainly show empty plots, this can be evidence of a strict motion rejection parameter, and optimization should be done.

Graphic Jump Location
Fig. 9
F9 :

(a) Block overlay layout shows number of free-noise epochs for each channel. (b) From these channels, several channels failed to extract any epochs due to heavy noise, as occurred in channel 1. (c) Channel 11 rejected four noisy trials while (d) there was no detected noise trial in channel 13.

Regarding the data-type conversion, POTATo will process continuous data following the listed Recipe order before the Blocking function and oppositely carry out block data processing after the Blocking function is applied. As described above, POTATo functions and layouts have specific data-type requirements. For example, Motion Check is merely functional for continuous data, while Band Filter works for both continuous and block data. Therefore, this may be noticeable when a user intends to modify the order of the Motion Check function. The Motion Check function cannot be listed after the Blocking function.

Baseline fitting

The resting (pre- and poststimulus) interval is assumed to show zero amplitude, yet an unstable baseline can occur due to systemic intervention delaying the baseline recovery.76,77 Considering undesirable matter, the Baseline correction function is applied to befit the data with adjustable fitting methods (e.g., averaged baseline, linear or polynomial regression) and subtract them into approximate zero-level. We set 1 s before stimulus onset and 5 s before epoch termination to be linearly fitted. This setting was chosen to include the effect of activation after stimulus such as second peak occurrence that was observed before in EDA. Therefore, to activate this function, we input 1-degree fitting and excluded the interval (stimulus onset 11  s after stimulus completion) in the parameter input-boxes. The baseline adjustment can be done not only in the block data but also in the continuous data using another function called Baseline fitting with the same principle; however, continuous fitting will result in a disadvantage in handling nonuniform baseline deviation, as depicted in Fig. 10.

Graphic Jump Location
Fig. 10
F10 :

Line plot of ΔCO2Hb signal in channel 22. (a) From continuous signal, baseline gradually decreased as well as performing task; however, sudden jump in baseline occurred and lingered. (b) Baseline in epoch data was corrected, suggesting approximate zero baseline.

After filtering and correcting the baseline, the cortical map layout and channel-wise plot are provided (plotting scroll-box → Plots and Map → Draw button), as shown in Fig. 11 (from another representative subject), to understand the distribution of ΔCO2Hb by selecting one signal type [Fig. 11(A)]. There are channels showing no expression (i.e., blank) and providing information that heavy noise had occurred and none of the epochs could be extracted from the filtering process. The interpreted map contouring the gap among channels can also be illustrated in this layout GUI [Mode scroll-box → Interpret; Fig. 11(B)]. If some channels are unavailable due to noise, the interpreted map connected to those adjoining channels will be impaired. Another feature, the time scroll button [Fig. 11(C)], is provided to display the time and assists the user in monitoring the ΔCHb alteration within the averaged epoch. According to the cortical map, the ΔCO2Hb signal noticeably enhanced in channels 12, 13, and 22 indicated by the red-colored gradient at 130 sampling points (i.e., 13 s, sampling rate 10 Hz) of the epoch.

Graphic Jump Location
Fig. 11
F11 :

Cortical map of 47 channels on prefrontal cortex area. Empty channels are caused by excessive noise. Map displays distribution of ΔCO2Hb at 13 s (sampling rate 10 Hz) after starting epoch or 12 s after stimulus onset time. Choosing one (A) signal type (B) and map type, visual observation over time course is done by (C) scrolling timepiece.

Statistical test of intertrials within subject

To confirm our hypothesis from the cortical map (Fig. 11) related to the enhanced ΔCO2Hb in channels 12, 13, and 22, we can average the ΔCHb amplitude and statistically verify the changes during the specified interval using the t-test function (All Filter scroll-box → t test → Add button). Figure 12(a) depicts the function GUI that enables a user to set both the tested baseline (Period 1) and stimulus-response interval (Period 2). In these WM tasks, periods 1 and 2 are determined as 10  s and 12–15 s (i.e., 0 s is onset stimulus time), respectively. POTATo will analyze the significance of the averaged ΔCHb in period 2 compared with period 1 with the data number depending on the total free-noise epochs for each channel. Because of the variability in peak time for each epoch,44 we can adjust additional interval time to find any peaks outside period 1 or 2 by clicking the Use Peak Search check-box. Furthermore, the independence level (α; 0.05) and result parameters can be set. After finishing with those parameters, OK button is clicked and we can examine the results in the cortical map layout (plotting scroll-box → Map for hdata.Results → Draw button). Figure 12(b) shows the results of the channel-wise statistical test with multiple statistical values (e.g., t-, p-, and h-values) for each data-kind (ΔCO2Hb, ΔCHHb, and ΔCHb-total). By selecting t-test Oxy.h [Fig. 12(A)], red channels denote p-values lower than α with positive t-values (enhanced ΔCO2Hb in period 2), while blue channel indicates p-value lower than α and negative t-value (enhanced ΔCO2Hb in period 1). The results confirm our hypothesis with other significantly enhanced ΔCO2Hb channels located adjacent to channels 12, 13, and 22 and the brain is seemingly activated in the right PFC. Nevertheless, this sole result is insufficient to determine ROIs and is also from only one subject. Hence, the grouped data are statistically analyzed to describe and estimate ROIs, as explained in Sec. 3.2.

Graphic Jump Location
Fig. 12
F12 :

(a) Parameter setting of intertrial t-test function within subject and (b) statistical results layout with selectable displayed parameters (A).

Additional features of preprocess graphical user interface

Before starting to create group analysis and statistically analyze the grouped data, we need to organize the Recipe for each data. If the dataset consists of a relatively large sample size, individually selecting functions and setting the parameters will waste much time. Therefore, POTATo provides the Recipe copy feature. By dragging several data, POTATo can detect the embedded Recipe and notify a user if there are any differences among Recipes, as shown in Fig. 13. This allows the user to choose one of the Recipes [Fig. 13(A)] and apply it to all selected data [Fig. 13(B)]. Additionally, the user can edit the function status and parameters in that Recipe [Fig. 13(C)]. In the Summary Statistical Analysis, the Recipe is assumed to be uniform or at least have the same data type; either continuous or block data. Now that we have already arranged the Recipe for all data, we can observe the averaged signal waveform by selecting several data and plotting them in Plot Average [Fig. 13(D)], as shown in Fig. 5(B), for spatial WM studies.

Graphic Jump Location
Fig. 13
F13 :

Different Recipe controller appears when several data are selected and different applied Recipes are detected. (A) Listed Recipes can be chosen to be applied for all (B) selected data. (C) This helps user to modify functions over subjects. After completing Recipe for each datum, (D) Plot Average provides average layout from all selected data.

For transforming the probe coordinate into the estimated Montreal Neurological Institute (MNI) space, a spatial analysis tool was developed according to the probabilistic registration7884 and is readily downloadable into the POTATo package.41 To start the spatial analysis, the Multifunctional button [Fig. 4(B)] is used by initially switching the button function to Spatial-Analysis (right clicking) and clicking that button after selecting the data. Once clicked, the Spatial Analysis DataViewer GUI appears. If the MNI panel is empty, the MNI template for a particular number of channels (e.g., 47 channels) can be loaded (P38\Benri Button\MNI\template\MNI_47ch.mat). By re-clicking the MNI panel, the three-dimensional (3-D) layout of the spatial registration is then provided together with the MNI coordinate and automatic anatomical labeling (AAL; i.e., Broadmann area estimation) for each channel. Instead of using Map for hdata.Results as the layout tool, MNI 3D hdata.Results is selected from the plotting scroll-box to visualize the statistical results in a 3-D registered-channel map.

In addition, the spatial analysis tool, POTATo provides the Project Repair feature in the Multifunctional button. As described in the POTATo concept section, POTATo enables a user to delete the Recipe via the POTATo Edit menu but not the POTATo Raw data. If the Recipe deletion mistakenly happened, the user can track down the missing Recipe by executing this Project Repair function. This will notify the user whether there are any existing Raw data without Recipes. Furthermore, the user can preferably determine to either create the Recipe or permanently remove the Raw data from that Project. Regarding the plug-in installment, the user can create the new plug-in from menu item (Tools → Plugin Wizard) by defining the plug-in type (e.g., preprocessing, summary statistical analysis, layout functions), used data-type (continuous and/or block data), and command code. This plug-in will be automatically installed in the POTATo package. A further explanation about plug-in creation is provided in the POTATo user’s guide included in the downloaded software package.

Data Grouping and Statistical Analysis

Before starting the experiment, we established the content and main objective as well as related hypotheses or expected results. The experiment was designed to accomplish the set goal and was conducted following this design. How to analyze and prove the initial hypotheses became major issues nonetheless. Statistical methods for resolving the getaway assist in the evaluation and standardization of experimental results and conclusion. To confirm the initial hypotheses and draw conclusions from a dataset, the computational statistical function can be performed. Preprocessed data are gathered, creating a summary of group analysis depending on the selected analysis approach in Summary Statistic Computation. Next, a statistical test (e.g., t-test) is carried out with conformable independent/target variables and paired/unpaired t-test. In the average sample-based approach, the average value of concentration change in the suspected activation interval is statistically analyzed instead of the estimated parameter for the predetermined explanatory HRF model (β) and fitting correctness (σ2 residual variance). Therefore, the average sample-based approach uses a single averaged value for each channel, as done in several previous studies,8588 while the model-based approach will take into account the data-model correlation and approximation (i.e., targeting minimized residual error through ordinary least square/OLS).

In a group analysis, the defective data may substantially affect the inferences. Therefore, subject screening is necessary as an initial step. For data validation, the channel is assumed to possess at least six free-noise epochs of eight conducted trials (75%) and each subject should retain at least 35 favorable channels. This was chosen to maximize utilization of dataset by 80%. By applying the similar Motion Check and Blocking recipes, the noisy epochs for each channel per individual subject can be identified. We recommend the Export WS button to extend usage outside the POTATo platform. By selecting the data with the same Recipe and clicking the Export WS button, the processed signals are computerized and transferred as workspace variables called data and hdata. The hdata facilitates easy retrieval because of hierarchical organization of specified information in a structure array including the marked noise called a flag (flag size=1 time point per channel), while the data variable stores the time-course signal of either continuous or block data depending on the Recipe. An epoch possessing even one flag is considered a defective trial and tends to be eliminated. Exporting more than a datum, POTATo will create structure arrays for each datum and store underneath the hdata cell array. By executing the code in the MATLAB® command window (Appendix 2), we can easily identify the number of defective epochs for each channel.

According to the predetermined criterion for 75% of valid epochs, the channels affected by noise in more than two epochs are categorized as defective. Therefore, the number of defective channels per subject is recognized. Figure 14 shows the relation between number of subjects and defective channels when performing a spatial WM task. For analyzing the number of subjects with controlled noise, we set utilization to at least 80% of the total number of subjects or 124 subjects. Based on noise data distribution, the subjects should have sustained more than 35 channels for the spatial WM task while 34 channels for the verbal WM task or around 74.5% and 72.3% of all channels, respectively. Therefore, we fixed the subject exclusion by 25 and 26 subjects for spatial and verbal WM tasks, respectively, and continued the analysis process in the Summary Statistics Computation and Statistical Test platforms. We now explain POTATo operation regarding two analysis approaches.

Graphic Jump Location
Fig. 14
F14 :

Graph of subject accumulation depending on channel validity (75% free-noise epoch per channel). For examining 80% of all subjects (128, orange line), up to 12 defective channels per subject is acceptable.

Average sample-based analysis

Prior to statistical analysis, the activation waveform should be observed to recognize the trend of ΔCO2Hb and ΔCHHb. Compared with introductory observation of averaged subject data [Figs. 5(A)], the waveform after sequent preprocess [Fig. 5(B)] mainly displayed a similar two-peaks pattern, indicating two types of activation. This is predicted as (1) encoding activation for memorizing the task and (2) recalling activation for retrieving the recent information. This assumption is related to the time course when the first peak appears during the delay of a given stimulus and the second peak appears after the memory test. This similar pattern might be interpreted as a positive sign of a tolerable preprocess without overdone correction. In addition, the preprocess step suggests a great advantage in heavy noise reduction. Consequently, the brain activity from noise-affected channels along the upper PFC can be extracted. Once again, from this waveform, the inference of the ROIs is still unclear.

Moving on to statistical analysis, the subjects from Data List are selected by excluding those with high noise-affected channels and accessing the Summary Statistical Computation platform by clicking the Summary Statistical Computation button located beside the Preprocess button [Fig. 4(D)]. Figure 15 illustrates the Summary Statistical Computation interface. To begin, a user is asked to name this group analysis in the Summary Statistical List input-box [e.g., Spatial/Verbal Task, Fig. 15(A)] and select one of the provided summary functions in the Function-Name scroll-box [Fig. 15(B)]. For the average sample-based approach, the Average function is selected. This function solely works for block data; the Blocking function must be included in the preprocessing Recipe. If some data are still processed in the continuous data, POTATo can identify them and the epoching feature will appear. In the Average function, there are three following features: subgroup analysis, baseline correction, and averaged period. First, the data characteristics listed in Key Search [Fig. 4(C)] can be used to divide a group analysis into smaller subgroups. POTATo also provides a more elaborate subgrouping that combines several characteristics in the Statistical Test GUI, the next interface. Second, if the Baseline correction function has not been incorporated into the preprocessing Recipe, POTATo allows the user to deliberately set this function parameter. Third, the activation interval evoked by the stimulus task together with the Peak Search feature can be adjusted similar to the t-test function. The initial waveform information addresses the concerned activation interval. In this analysis, the activation interval was selected on the basis of one-fourth wavelength; thus, the intervals of 6 to 8 s and 12 to 15 s after stimulus onset were determined as activated for the first and second peaks, respectively. Because the Motion Check function is also included in the processing Recipe, the Apply Flags check-box is clicked to use the information of defective epochs (flag) and leave out those epochs in further analysis. Besides the Function-Name box, there are additional buttons such as Save, Load, and Draw. The Save and Load buttons work for storing the information of adjusted parameters in a .mat file and launching the saved information to other group analyses, respectively. The Draw button also provides the layout function similar to the Plot Average button. This setting is then saved to proceed with the statistical analysis.

Graphic Jump Location
Fig. 15
F15 :

Summary Statistics Computation interface creates multiple group analyses with (A) specific name and (B) function. While editing former group analysis setting, data with modified Recipe is notified with (C) three asterisks and enables user to select either using previous Recipe (‘Use Summary Statistic Data Recipe’ command) or newly modified (D) Recipe (‘Use Analysis Data Recipe’ command).

As discussed earlier, the degree of amplitude signal can indicate the activated ROIs. To minimize the intravariability signal and ensure robust analysis, normalized amplitude is calculated by comparing the averaged stimulus and rest amplitudes with the following equation for each channel34Display Formula

where X is the mean-amplitude of either ΔCO2Hb or ΔCHHb from free-noise epochs, SD is standard deviation, and n is the number of free-noise epochs in channel-wise manner. In the Average function, the concerned interval is merely the evoked activation interval. Therefore, the t-test function is applied for each data subject with an averaged period according to the predetermined rest/stimulus interval to obtain this parameter magnitude. The normalized value is stored in a Results structure under hdata; thus, the group analysis is created in the Summary Statistics Computation platform by using another function, save hdata.Results.

To access the Statistical Test GUI, the user clicks the Statistical Test button located beside the Summary Statistics Computation button [Fig. 4(D)]. As illustrated in Fig. 16, t-test is performed by initially selecting the summarized group analysis [Summary Statistic Data List scroll-box → ‘the recent group analysis’ → Add button; Fig. 16(A)]. There are several statistical test functions such as t-test for average sample-based and model-based approaches [Fig. 16(B)]. The customized analysis parameters, such as independent and target variables, allow the user to deliberately explore their own data [Figs. 16(C) and 16(D)]. To create elaborate subgroups as described above, the Target Group feature is used allowing the user to sort several data characteristics (e.g., signal type, filename, age, gender, and so on) and form a target group [Add button; Fig. 16(E) → choose characteristic grouping → OK button]. To conduct two-sample analyses, the user has to create two target groups. For instance, to understand the differences between spatial and verbal WM tasks, the first target group is set to use the ΔCO2Hb signal type from spatial task, whereas the second target group uses the same ΔCO2Hb signal type and sources from verbal task. To analyze in a channel-wise manner for each WM task, channel, normalized amplitude (i.e., z-value), and ΔCHb signal type (e.g., ΔCO2Hb or ΔCHHb) are set as the independent variable, target data, and target group, respectively. Furthermore, the average target data feature is set beforehand [Fig. 16(F)] to compute either each trial or intertrial average as one datum. The Launch t-test button [Fig. 16(G)] is then clicked to display another analysis GUI enabling the user to modify the statistical hypothesis parameters, such as independence level and tail type, to compute several analysis outputs (i.e., t- and p-values), and to display the analysis output in a spatially registered-channel map as illustrated in Fig. 17. For further analysis in the MATLAB® window or other statistics software, the Text Out and Export WS buttons [Fig. 16(H)] can be used to export the summary of group analysis in .csv and .mat files.

Graphic Jump Location
Fig. 16
F16 :

(A) Statistical Test interface with several setting parameters such as selection of group analysis, (B) statistical analysis function, (C) independent variable, (D) target data (summary of group analysis), (E) target group, and (F) optional subaveraged target data. For executing analysis, (G) launch t-test button is clicked to perform analysis in different GUI. (H) Exporting tool is provided to analyze summary of group analysis outside POTATo package.

Graphic Jump Location
Fig. 17
F17 :

Significance results (one-sample t-test with right tail for ΔCO2Hb and left tail for ΔCHHb) of resting-activation period for (a, b, e, f) spatial and (c, d, g, h) verbal WM tasks with respect to (a, c, e, g) ΔCO2Hb and (b, d, f, h) ΔCHHb (a)–(d) signals of encoding and (e)–(h) recalling periods. p-values are shown in either negative or positive logarithmic scale for ΔCO2Hb and ΔCHHb results, respectively. Marked regions (c, d, g, and h) show probable left laterality in verbal task while recalling memory.

In general, the evoked neural activity is associated with the elevated energy demand through glucose catabolism coupling with increased neurometabolic and neurovascular responses.31,89 Consequently, the activated cortex expresses incremental cerebral blood flow (CBF) together with cerebral blood volume (CBV) delivering oxygen and nutrients to compensate for the neural-task activity. While a higher metabolism rate is promoted in parallel, the oxygen supply from CBF and CBV is greater; therefore, deoxygenated hemoglobin encounters concentration depletion. Nevertheless, the neurometabolic and neurovascular responses are more complex than the widely known theory especially at non-ROIs that are sometimes excluded in analysis. For example, there are some channels at the bilateral superior temporal gyrus, premotor cortex, and frontopolar areas (i.e., channels on left, right, and lower probe perimeters) where both increased ΔCO2Hb and ΔCHHb signals are shown in Fig. 5(B). Several hypotheses, such as influential physiological noise at skin blood flow with both increased ΔCO2Hb and ΔCHHb90 or spatial inhomogeneity of motion artifact (e.g., eye blinking or frowning may strongly influence the frontopolar compared with other regions) with enhanced positive correlation between ΔCO2Hb and ΔCHHb,25 have been introduced; however, further investigation is necessary to address this issue.

To evaluate specific ROIs over measured regions, the multiple hypotheses tests are involved and increase the familywise error rate because of incorrectly rejecting the null hypothesis of no evoked-stimulus activation (i.e., false discovery or type-1 error).9193 This multiplicity has been an important and inevitable issue in any neuroimaging techniques; thus, the Bonferroni correction (α=0.01/47) was adopted to avoid type-1 error and maintain the familywise error rate in analyzing these WM tasks. According to the results shown in Fig. 17, the ΔCO2Hb concentration apparently increased in wide areas regardless of task types (p<0.00022) presumably confirming the insufficiency of statistical correction (i.e., Bonferroni) or noise elimination. These increases again raise the question of particular ROI positions. Nonetheless, there were some channels showing distinguished significance degrees with extremely low p-values [red-colored channels in Figs. 17(a), 17(c), 17(e), and 17(g)]. Extensively, the significance of decreased ΔCHHb promoted ROI information as a hypothesis of opposite correlation with incremental ΔCO2Hb [blue-colored channels in Figs. 17(b), 17(d), 17(f), and 17(h)]. Simultaneously checking the p-distribution map for both task and signal types, channels 13, 16, 22, 26, 32, and 35 were more likely to behave as ROIs and those channel locations reflected to each other (channels 13–16, 22–26, and 32–35) through the vertical nasion axis. Expectedly, those channels had also been predicted from the initial visual observation (EDA).

According to the AAL of the spatial analysis tool, those channels are mainly estimated as the bilateral dorsolateral PFC (DLPFC; Broadmann area 9 and 46), pars triangularis Broca’s area [partly ventrolateral PFC (VLPFC), Broadmann area 45], and tiny section of the frontopolar PFC (Broadmann area 10), as illustrated in Fig. 18. These results confirmed the similar ROIs from previous fNIRS and fMRI studies related to activated-WM ROIs.9496 Despite the EDA and MD results of increased ΔCO2Hb in the broad frontopolar PFC area (Broadmann area 10), the hypotheses of systemic interference and inhomogeneous noise affecting particular regions might contaminate the actual activation signal. This was also confirmed by the significantly increased ΔCHHb, which went against the negative correlation with ΔCO2Hb at the activated regions. In contrast to the bilateral activation of the spatial stimulus, the verbal WM task differed in probable left-laterality relying on a greater agreement area between increased ΔCO2Hb (Bonferroni corrected) and decreased ΔCHHb (uncorrected) criteria in the left DLPFC [Figs. 17(c), 17(d), 17(g), and 17(h) with marked regions], especially in suggested recalling process. In spite of reported right- and left-laterality for spatial and verbal WM tasks,9799 respectively, the left-preference of a verbal task was solely inferred in these datasets.

Graphic Jump Location
Fig. 18
F18 :

Spatial estimation of activation area (yellow channels) with indicated particular Broadmann area (BA) at respective channels (Channel 13: BA9 98.95%, BA46 1.05%; Channel 16: BA9 90%, BA46 10%; Channel 22: BA45 37.18%, BA46 62.82%; Channel 26: BA45 68.6%, BA46 31.4%, Channel 32: BA10 16.28%, BA46 83.72%; and Channel 35: BA10 6.17%, BA46 93.83%).

Model-based approach using general linear model analysis

Other than the mean value, the evaluation of a continuous signal is performed by the GLM estimating the likelihood of the HRF model. As a successor of the GLM implementation in the fNIRS data, Minagawa-Kawai et al.100 established the GLM analysis plug-in for POTATo platform that works in four steps. First, the establishment of a group analysis and execution of the GLM using the boxcar function as multiple regressors in the Summary Statistics Computation interface (installed keio-lscp GLM function).101 Second, the finite impulse response (FIR) method is implemented on the boxcar regressors to determine the optimum time latency of the HRF model via the keio-lscp GLM test function in the Statistical Test interface.102 Third, the GLM is re-performed using the HRF model with obtained time latency in the keio-lscp GLM function (Summary Statistics Computation interface) and the β parameter is further estimated. Finally, the brain activity level indicated by the acquired β value is evaluated using the t-test (keio-lscp) function in the Statistical Test interface. While the free plug-in and user’s guide are available on their website (Japanese only),103 this tutorial also provides a brief explanation of plug-in incorporation into the POTATo interface.

To get started, this free plug-in (P3_keio-lscpGLM_v1.02.zip) is downloaded and the extracted file consists of three folders called PluginDir, Research1stFunctions, and Research2nd Functions. Subsequently, the contents of these folders are copied into the POTATo package (P3) with matching folder names. POTATo is re-launched and the user can operate this plug-in. After selecting the data from Data List and switching the functional interface to Summary Statistics Computation, the keio-lscp GLM is selected instead of the previous Average function [Fig. 15(B)]. There are several input parameters, such as epoching interval, GLM, error parameter regressor (e.g., sin and cos), and ΔCHb types, as illustrated in Fig. 19. Because this function also includes epoching [Fig. 19(A)], preprocessing is done on the continuous data and the Blocking function is removed from the Recipe. For the first step, Boxcar [Fig. 19(B)] is selected as multiple regressors with a modifiable regressor width, number, and asynchrony interval. We set those parameters to 1.0-s width, 20 boxes, and 1.0-s asynchrony interval of regressors. By combining Boxcar and sin-cos error regressors, the Check Design Matrix button [Fig. 19(C)] provides visualization of the design matrix for the concerned participant, ΔCHb types and channel, as depicted in Fig. 20(a), where the x- and y-axes denote the length of the sampling point and number of regressors, respectively. If there are any defective sampling points, the number of regressors will also increase to compensate for the detected noise. The setting for the ΔCO2Hb and ΔCHHb data is done separately.

Graphic Jump Location
Fig. 19
F19 :

keio-lscp GLM function in creating group analysis in Summary Statistic Computation interface. Epoching feature is (A) embedded with selectable model to design (B) adaptive design matrix. Channel-wise design matrix can be observed for each data-kind (ΔCO2Hb, ΔCHHb, and ΔCHb-total) and (C) each subject.

Graphic Jump Location
Fig. 20
F20 :

Established design matrix using (a) boxcar function and (b) optimized HRF-time latency model.

Second, the optimization of time latency fitted to the canonical HRF model is accomplished in the Statistical Test interface by applying the keio-lscp GLM test function, as illustrated in Fig. 21(A). In this GUI, the user is allowed to modify the variables (File and Channel scroll-boxes and stimulus marker, Fig. 21(B) to optimize the time latency. If the user would like to determine the optimum condition from the average of all subjects and channels, all subjects (File scroll-box) and channels (Channel scroll-box) are selected (CTRL+A) and both scroll-boxes are set as the Mean. For channel-wise analysis, the time latency for each channel is optimized by averaging the respective channels from all subjects. Therefore, after selecting all subjects and channels, the File method is set as the Mean, while the Channel method is defined as None. By clicking the Check FIR button [Fig. 21(C)], the plot layout according to the set parameters will appear. Fitting to the HRF model is done afterwards by loading the canonical model file (P38\PluginDir\keioGLM\keioGLM_HRF32s.mat), and the optimization results are examined by clicking the Check button before saving it in the .mat file [Fig. 21(D)].

Graphic Jump Location
Fig. 21
F21 :

(A) Statistical Test interface for optimizing time latency of HRF model using keio-lscp GLM test. (B) Optimization can be adjusted depending on set parameters such as intrasubject and channel-wise averages. (C) Check FIR button provides visualization of output signal before being fitted to HRF model. (D) Optimization result can be confirmed and needs to be saved for re-using GLM with optimized time latency model.

Figure 22 represents the approximation results of the fitted HRF model from the average channel-wise signals of all subjects during the spatial WM task. It should be noted that the exploratory waveform pattern is heretofore observed as the appearance of two peaks in a single time task trial. It is presumably disparate from the common HRF model with one peak, oscillating over, and returning back to the rest state. Even though manipulation can be done by decreasing the number of regressors to insensitively scan the waveform, it is unwise to do this for primarily satisfying the HRF-model likelihood. Together with the misleading discovery of the activated frontopolar PFC in the preliminary averaged-raw signal, these prove the importance of EDA-MD parallelism, with which the user should not visually over-interpret the qualitative results or aimlessly follow the known model and ignore the uniqueness of the data.

Graphic Jump Location
Fig. 22
F22 :

(a) Representative model comparison among HRF, actual waveform of spatial WM task (FIR), and corrected-delay approximate HRF models. (b) Optimization of time latency shows delay time of 4.4 s. (c) HRF model with optimized time latency mode.

To re-perform the GLM using optimized time latency, we return to the Summary Statistics Computation and edit the current group analysis (e.g., Spatial/Verbal Task). Instead of using the Boxcar function, Optimal Model from File made in Stat is chosen from Use Model scroll-box [Fig. 19(B)] and the saved optimization result (.mat) is loaded. Similar to the first step, the design matrix for each channel and participant can be confirmed (Check Design Matrix button), resulting in different matrices consisted of the HRF model, fluctuation noise, and constant element regressors [Fig. 20(b)] since the boxcar function and optimum condition are performed. Clicking the Apply changes button revises the summary of group analysis made from the Boxcar function. To privately analyze the β parameters outside the POTATo interface, the Text Out button can be used to import the estimated parameters in .csv format. The fourth step is then performed in the Statistical Test interface to statistically verify whether the β parameter is significantly greater than 0. By selecting the t-test (keio-lscp) function [see. Fig. 21(A)], the user is allowed to customize the independent/target group variables and sample analysis type (e.g., one- and two-sample test). Similar to the general t-test function, the Perform t-test button provides another analysis GUI to set the hypothesis parameters, such as independence level (α=0.01/47) and tail type, and to execute the statistical analysis. The 2-D statistical result map can be transferred into a 3-D spatially registered map by exporting the result variables (e.g., p- and t-values) and assigning them to the spatial analysis tool, as described above.

Even though the GLM method was restricted by the inflexible HRF model that overlooks the variability of task-response model, the results of t-test with Bonferroni correction (Fig. 23) showed similarity with the previous average sample-based approach. This might have been caused by the incorporated peak interval within 20-s regression. However, the particular effects of each peak representing the different activations of encoding and recalling could not be examined separately. Analogous to ROIs determination, these results were analyzed based on the primary hypothesis of ΔCO2Hb and ΔCHHb correlation. By concurrently analyzing both averaged sample and GLM results in not only spatial but also verbal WM tasks, there seemed to be significant overlapping and neighboring channels in both ΔCO2Hb and ΔCHHb. Further, the ROIs were obtained at nearly identical channels with greater left PFC predilection for both WM tasks (Fig. 23 with marked regions) verifying the analysis reproducibility. From these results, the formerly selected ROIs were reasonable and also confirmed the preceding results of WM-related ROIs at DLPFC.

Graphic Jump Location
Fig. 23
F23 :

Activation distribution map (one-sample t-test, one-tail) of HRF-likelihood model for spatial (a, b) and verbal (c, d) WM tasks. Both ΔCO2Hb (a, c) and ΔCHHb (b, d) signals conclusively present close tendency with activated region in averaged-amplitude assessment. Marked regions show probable left-laterality for spatial and verbal WM tasks (p-value in negative logarithmic for ΔCO2Hb and oppositely for ΔCHHb).

Additional feature of Summary Statistics Computation GUI

Due to trial-and-error attempts in analyzing a dataset, we frequently switch the Preprocess—Summary Statistics Computation—Statistical Test interfaces. While Recipes are modified which affect the processed signals, the summary of group analysis remains the same. Before performing Statistical Test, the summary of group analysis should be edited and updated. To inform the user, the Data List of the Summary Statistics Computation interface is automatically marked with three asterisks (***) if preprocessing Recipes are modified, as illustrated in Fig. 15(C). By selecting all data and right-clicking the Data List, there is another notification [Fig. 15(D)] requesting the user to either accept the recent Recipe modification or reject the modification. In case of accepting modification, the user can update the summary of group analysis by applying the new Recipe and clicking the Apply Changes button [Fig. 15(E)]. In contrast, rejecting modification will alter the preprocessing Recipe into the former Recipe saved in the summary of group analysis.

POTATo as an IDEA platform provides the flexibility of preprocessing, analyzing, and visualizing data. POTATo is expected to assist a user in optimizing functional brain activity interpretation by supporting the parallelism between EDA and MD analysis. Not only constantly examining data with conventional approaches, it is also possible to establish and develop analysis methods in POTATo. The integration of sequential analysis steps in a platform is expected to enhance data analysis comprehension and promote broader fNIRS application. By using WM tasks as analysis examples, we discussed POTATo’s potential in performing analysis according to average sample-based and model-based approaches. Initiated with EDA of raw data and confirmed through MD evaluation of processed data, these results validated previous studies on activated DLPFC while performing various WM tasks (delay-time, N-back tasks). This also confirms the analysis practicability of the POTATo platform.

For future fNIRS studies, there are issues that should be addressed, mainly regarding noise, (1) the noise source and identification,13,104106 (2) perplexity of noise response by either rejecting the trial or correcting the noise-affected signal with optimized method,107109 and (3) the diverse distribution of noise in spatial and temporal frames.4446,110112 Overcoming these issues will become the current trend in fNIRS studies. Not only noise elimination in data analysis but also improved noise prevention in hardware (e.g., multidistance optodes)113116 and experimental design should be done simultaneously. As well as helping users develop an analysis method to solve those issues, we are also continuously improving POTATo’s features to provide a more convenient and resourceful analysis domain.

POTATo is able to import both intensity and ΔCHb data allowing users to provide one of those. Tables 1 and 2 are the examples of data structure in .xls file format. The necessity of wavelength information depends on the availability of raw data type which the ΔCHb data does not require the wavelength information. Besides the wavelength, other information, such as file name, ID, subject name, age, gender, measurement date, sampling rate, stimulus mode, measurement data, and stimulus-related marker, are compulsory as the minimum requirements of POTATo import function. Table 3 presents the details of each component in the data structure.

Table Grahic Jump Location
Table 1Example of intensity data.