3.1.1.

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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}/\mathrm{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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ signal (channel 11, $3\times 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ amplitude spectrum (upper plot) that might be related to the $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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.

Fig. 6

(a) Spectrum Check layout in representative $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ signal from channel 11 while performing spatial WM task. Upper plot shows $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ amplitude spectrum versus frequency and lower plot illustrates $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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,⁷³ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$, $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, and $\mathrm{\Delta }{C}_{\mathrm{Hb}\text{-}\text{total}}$, respectively.

Fig. 7

Motion Check interface conveniently confirms noise occurring at particular time points for each channel and data-kind (1 for $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$, 2 for $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, and 3 for $\mathrm{\Delta }{C}_{\mathrm{Hb}\text{-}\text{total}}$). Impact on channel can be immediately monitored as parameter justification. Detected $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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\mathrm{mM}·\mathrm{mm}$, kind 1, and 2 points in the threshold, data-kind, and interval point boxes, respectively, POTATo will look into the ButterWorth-filtered $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ signal for each channel, detect whether there are any $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ deviations over the $0.1\text{-}\mathrm{mM}·\mathrm{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\mathrm{mM}·\mathrm{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.

3.1.2.

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,¹⁸ 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., $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$, $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, and $\mathrm{\Delta }{C}_{\mathrm{Hb}\text{-}\text{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 $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ 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.

Fig. 8

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 $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ (C) and (D) reconstruction plots for Raw data (A) and (C) and filtered signals (B) and (D).

3.1.3.

Task-related signal

POTATo was designed by adopting the block design analysis^74,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.

Fig. 9

(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.

3.1.4.

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\mathrm{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.

Fig. 10

Line plot of $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ alteration within the averaged epoch. According to the cortical map, the $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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.

Fig. 11

Cortical map of 47 channels on prefrontal cortex area. Empty channels are caused by excessive noise. Map displays distribution of $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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.

3.1.5.

Statistical test of intertrials within subject

To confirm our hypothesis from the cortical map (Fig. 11) related to the enhanced $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ in channels 12, 13, and 22, we can average the $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ 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 $-1-0\mathrm{s}$ and 12–15 s (i.e., 0 s is onset stimulus time), respectively. POTATo will analyze the significance of the averaged $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ 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,⁴⁴ 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 ($\alpha$; 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 ($\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$, $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, and $\mathrm{\Delta }{C}_{\mathrm{Hb}\text{-}\text{total}}$). By selecting $t$-test Oxy.h [Fig. 12(A)], red channels denote $p$-values lower than $\alpha$ with positive $t$-values (enhanced $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ in period 2), while blue channel indicates $p$-value lower than $\alpha$ and negative $t$-value (enhanced $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ in period 1). The results confirm our hypothesis with other significantly enhanced $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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.

Fig. 12

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

3.1.6.

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.

Fig. 13

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 registration^78–84 and is readily downloadable into the POTATo package.⁴¹ 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.

3.2.1.

Average sample-based analysis

Prior to statistical analysis, the activation waveform should be observed to recognize the trend of $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$. 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.

Fig. 15

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 channel³⁴

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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ signal type from spatial task, whereas the second target group uses the same $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ signal type (e.g., $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ or $\mathrm{\Delta }{C}_{\mathrm{HHb}}$) 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.

Fig. 16

(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.

Fig. 17

Significance results (one-sample $t$-test with right tail for $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and left tail for $\mathrm{\Delta }{C}_{\mathrm{HHb}}$) 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) $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and (b, d, f, h) $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ (a)–(d) signals of encoding and (e)–(h) recalling periods. $p$-values are shown in either negative or positive logarithmic scale for $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ signals are shown in Fig. 5(B). Several hypotheses, such as influential physiological noise at skin blood flow with both increased $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$⁹⁰ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$,²⁵ 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).^91–93 This multiplicity has been an important and inevitable issue in any neuroimaging techniques; thus, the Bonferroni correction ($\alpha =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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{\mathrm{C}}_{\mathrm{HHb}}$ promoted ROI information as a hypothesis of opposite correlation with incremental $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ [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.^94–96 Despite the EDA and MD results of increased $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, which went against the negative correlation with $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ (Bonferroni corrected) and decreased $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ (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,^97–99 respectively, the left-preference of a verbal task was solely inferred in these datasets.

Fig. 18

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%).

3.2.2.

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.¹⁰⁰ 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).¹⁰¹ 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.¹⁰² 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 $\beta$ parameter is further estimated. Finally, the brain activity level indicated by the acquired $\beta$ 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),¹⁰³ 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 $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ 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, $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ data is done separately.

Fig. 19

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 ($\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$, $\mathrm{\Delta }{C}_{\mathrm{HHb}}$, and $\mathrm{\Delta }{C}_{\mathrm{Hb}\text{-}\text{total}}$) and (C) each subject.

Fig. 20

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)].

Fig. 21

(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.

Fig. 22

(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 $\beta$ 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 $\beta$ 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 ($\alpha =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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ 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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$. 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.

Fig. 23

Activation distribution map (one-sample $t$-test, one-tail) of HRF-likelihood model for spatial (a, b) and verbal (c, d) WM tasks. Both $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ (a, c) and $\mathrm{\Delta }{C}_{\mathrm{HHb}}$ (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 $\mathrm{\Delta }{C}_{{\mathrm{O}}_2\mathrm{Hb}}$ and oppositely for $\mathrm{\Delta }{C}_{\mathrm{HHb}}$).

3.2.3.

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.