OptoNet II: An Advanced MATLAB-Based Toolbox for Functional Cortical Connectivity Analysis With Surrogate Tests Using fNIRS

Cortical connectivity analysis is a widely used method for understanding the causes of neurological disorders and related brain mechanisms. Although there exist numerous activity analysis toolboxes for functional near-infrared spectroscopy (fNIRS), there are only a few cortical connectivity analysis toolboxes. In 2019, we released a MATALB toolbox named OptoNet, which has helped researchers to analyze brain networks using fNIRS. In this study, we developed an advanced MATLAB toolbox, named OptoNet II, to add new features that overcome the shortcomings of OptoNet. With these new features, OptoNet II can efficiently analyze cortical connectivity according to brain region using any fNIRS channel sets and can present the results of two connectivity analyses with auto-thresholding based on surrogate tests. To evaluate the efficacy of the new functions, the finger-tapping task experiment was carried out before and after transcranial direct current stimulation (tDCS) in the primary motor area. OptoNet II can efficiently show the effects of tDCS on functional brain region connectivity, which has been difficult to confirm by conventional methods. In this article, we propose the OptoNet II as a useful and efficient toolbox for researchers who want to perform cortical connectivity analysis using fNIRS.


I. INTRODUCTION
Functional neural connectivity plays an important role in advanced cognitive processes, and has thus drawn increasing attention from researchers over the past decades [1]- [6]. Various measures, such as phase synchronization index (PSI) [7], mutual information [8], partial directed coherence [9], frequency ratio [10], and mean phase coherence [9], have been applied to quantify the functional connectivity of different brain units using multichannel neural signals, including electroencephalography (EEG), functional magnetic resonance imaging (fMRI), and functional near-infrared spectroscopy The associate editor coordinating the review of this manuscript and approving it for publication was Poki Chen .
(fNIRS) [5]. These functional connectivity measurements have shown similar global functional connectivity patterns across studies, with some differences between study results for particular cortical regions, and have also shown agreement with regards to quantifying the level of synchronization [2].
fNIRS is a noninvasive method used to measure hemodynamic brain signals based on the absorption of near-infrared light, with wavelengths in the range of 650 nm to 950 nm, transmitted through the intact skull [11]. fNIRS monitors variations in regional cerebral blood flow and estimates hemodynamic signals, which are highly correlated to the blood-oxygenation-level-dependent (BOLD) signal outputs in fMRI [12]. The important advantages of fNIRS are its low cost, portability, and the potential to extend research to a VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ wider range of environments than many other neuro-imaging systems, such as fMRI, EEG, and magnetoencephalography (MEG) [13]. General linear modeling (GLM) is one of the most extensively used models to represent data in a linear combination form, and constitutes a standard method for analyzing fMRI data. Indeed, many statistical analysis toolboxes have been developed for fNIRS based on the GLM. However, GLM-based analysis methods often fail to adequately analyze brain functions because of artifacts in the fNIRS measurements. The artifacts exist for various reasons, such as subject movements, blood pressure variations, and instrumental instabilities [14]- [19]. Recently, various connection and causality estimation methods for functional brain network analysis have been developed, and have demonstrated their utility in cognitive neuroscience and neurological clinical studies [20]. Many brain network estimation methods have been applied to numerous functional neuroimaging modalities, such as EEG, local field potential, intracranial EEG, MEG, fMRI, and fNIRS.
The most popular tools currently available for fNIRS analysis are hypergeometric optimization of Motif enrichment (HOMER) from Harvard Medical School [21] and near-infrared spectroscopy statistical parametric mapping (NIRS-SPM) from the Korea Advanced Institute of Science and Technology (KAIST) [15]. HOMER (available at http://www.nmr.mgh.harvard.edu/PMI/) calculates individual hemodynamic responses using ordinary least-squares linear deconvolution, while NIRS-SPM (available at http://bisp.kaist.ac.kr/NIRS-SPM) applies the SPM method, which refers to the construction and assessment of spatially extended statistical processes used to test hypotheses about functional imaging data [3], [15]. However, these analysis tools cannot estimate functional brain connectivity, and they are difficult for new users to learn. Therefore, prior to the release of OptoNet, there was no software program for the analysis of functional brain networks using fNIRS available for free and which could be used by unskilled users [3]. In 2019, our research team released a Windows-based graphical user interface (GUI) MATLAB toolbox named OptoNet with the aim of providing unexperienced users with an easy way to analyze cortical networks based on 3-dimensional (3D) finite element analysis (FEA) [22], [23]. However, OptoNet can only analyze cortical networks for each channel of any given fNIRS system; thus, if the system has many channels, the results of the cortical network analysis can be too complex to check certain networks. Furthermore, OptoNet uses a manual threshold setting, which can only show those networks over the user-determined threshold value. This manual threshold setting can cause the objectivity and reliability of the results to be poor.
In this article, we introduce an advanced version of OptoNet, named OptoNet II, which is also a Windows-based MATLAB toolbox. The primary differences between OptoNet II and OptoNet are the cortical connectivity analysis according to the functional brain region, the auto threshold method that uses a surrogate test, and the representation function of the connectivity difference.
In contrast to the previous version of OptoNet, users of OptoNet II can freely set all fNIRS channels to fit the functional areas of the brain, and therefore analyze cortical connectivity according to functional region. Additionally, OptoNet II features an auto-threshold function based on surrogate tests, which can provide significance tests if the original signal has its property randomized among the surrogate data [10]. Detailed descriptions and examples of OptoNet II are provided in the following sections. We tested OptoNet II using 64-bit Windows 10 installed on Intel i5 and i7 personal computer systems.

II. METHODS
The procedures of OptoNet II can be roughly divided into the following steps: signal processing, brain region model setting, and connectivity analysis. In the fNIRS signal processing step, loading fNIRS data, selecting the fNIRS epoch, and setting the analysis duration are performed. The connectivity processing step consists of loading the standard head model, setting the brain regions, and executing connectivity analysis. The GUI of the toolbox is optimized for the latest version of the MATLAB (ver. R2020a) (MathWorks, Natick, MA, USA) and GeForce (Nvidia, Santa Clara, CA, USA) graphics card series. Figure 1 shows the overall sequential steps and GUI of OptoNet II. Next, detailed descriptions for each step will be presented.

A. SIGNAL PROCESSING STEP
The signal processing step is performed in the signal processing panel of OptoNet II. First, the measured fNIRS data can be entered into the ''NIRS Data Load'' section indicated in Figure 2 A. OptoNet II supports the.NIR file format, which can be easily converted from the.mat file format of MATLAB. OptoNet II also provides an NIRX converter, as shown in Figure

B. BRAIN REGION MODEL SETTING STEP
As seen in Figure 3, the connectivity processing step is performed in the connectivity processing panel. Users can select between the standard and the individual head and cortex models in the ''Load Head and Cortex Model'' section shown in Figure 3 A. Figure 3 B shows the ''Region Analysis'' section. If the on/off checkbox is checked, the region analysis mode is activated, while if it remains unchecked, OptoNet II is essentially in fNIRS channel-based analysis mode, similar to the original OptoNet. Users can set up brain region modeling according to fNIRS channels by inputting information  Once entered, the number of fNIRS channels for each functional region can be saved and used again later. The saved brain region model can be loaded by using the button in Figure 3 B, and then the brain regions that will be analyzed are represented in Figure 3 F. If the user wants to make a custom brain region model that has a different number of functional brain regions from the standard region, the user decides the positions of the region maker on the head model using the position setter in Figure 3 D. Then, custom fNIRS channels are entered as many as the number of the custom regions using the custom region setter that repeats pop-up until pressing the save button in Figure 3 E. The fNIRS channels that are set to the same functional brain region are grouped into that region, and connectivity analysis is performed according to the set brain region model.

C. CONNECTIVITY ANALYSIS STEP
One of the biggest problems that exists when estimating brain connectivity is determining whether corresponding brain region pairs are significantly connected or not with high confidence. The previous version of OptoNet only offered the manual threshold setting, which could only represent connections over the manually set threshold; therefore, it could not determine whether there was significant connectivity between regions. In order to overcome this problem, the surrogate data can provide significance tests on whether the original signal has the property randomized among the surrogate data [24], [25]. Therefore, an artificial surrogate test was adopted in this version of OptoNet II.
Surrogate methods usually produce artificial data by randomizing the property to be tested while mimicking other properties (e.g., the spectra) of the original signal as much as possible; the more properties that are preserved, the stricter the generated surrogate data [26]. The surrogate time series has the same rank sequence as the Gaussian time series, but the samples of the surrogate time series all come from the measured fNIRS signal. This means that the surrogate series is a rank-shuffled version of the fNIRS signal, with the temporal structure destroyed but the distribution, mean, and variance all preserved [26]. In other words, the surrogate test randomly and independently shuffles fNIRS signals from each cortical source to create a surrogate data set. In order to shuffle fNIRS signals, the phase and amplitude were randomized in the Fourier transformed fNIRS signals to preserve the power spectra feature of measured fNIRS signals with less distortion. Network analysis was estimated from the surrogate data set by performing the same calculating process for 10,000 or more iterations (10,000 times is a default value in the software). In this way, the empirical distribution for a given network estimator can be created. From the distributions, null hypotheses are provided: ''the independently and randomly generated surrogate data set has no interactions between generated fNIRS signals''. From these distributions, significant levels for estimated networks were determined from each empirical distribution [2], [27]. The surrogate tests used the auto-threshold method that was added in this new version of OptoNet II, and which can be activated using the actions shown in Figure 3 G. To select significant networks, the toolbox follows a rank-order test [28]. First, a residual probability α of false rejection is selected corresponding to a level of significance ((1 − α) × 100%) [25]. Then, for a one-tailed test, which can look for small prediction errors, surrogate sequences (M = K/α − 1) are generated, where K is a positive integer. For a two-tailed test, which can go both ways for time asymmetry, the surrogates (M = 2K/α − 1) are generated, resulting in a probability α that the data gives either one of the K smallest or largest values [25]. Using more surrogates can increase the discrimination performance [25], and 10,000 surrogate samples were used in OptoNet II. If a network has a higher surrogate than the level of significance, it will be represented as a significant network.
The method used for the connectivity analysis can be selected from one of following: correlation [29], [30], coherence [31], frequency ratio [32], phase locking value (PLV) [33], and is executed as shown in Figure 3 G. After analysis, only significant connectivity is represented in Figure 3 F by auto-thresholding based on the surrogate tests. The resulting values for all estimated connectivities, including insignificant connectivities, are saved in a work folder as a.mat file. If the user wishes to compare two connectivity result files (e.g., rest state and task state), the differences in connectivity can be estimated by using the button shown in Figure 3 H. The difference in the connectivities is calculated through a simple subtraction from the two estimated connectivity results, and then the significant connectivities of the difference are selected using the surrogate test again. After that program is executed, the two results files are continuously loaded and the differences in connectivity between the two files are calculated. Then, the significant differences in connectivity are graphically represented on the head model shown in Figure 3 F.

III. EXPERIMENTAL RESULTS
We performed an experiment that highlights the improved features of OptoNet II over the original version. The experimental task paradigm consisted of the finger tapping task using the left hand. In this experiment, finger tapping is performed sequentially, and participants must accurately tap their finger when prompted by a sign displayed on the monitor screen shown in the right picture in Figure 4 A. Twenty-three adult volunteers over 19 years old who were right-handed were recruited. All subjects were healthy without a history of brain injury, neurological, or psychiatric diseases. Study protocols were approved by the Samsung Medical Center (SMC) Institutional Review Board (IRB), and all subjects gave written informed consent prior to each experiment.
A. fNIRS SYSTEM fNIRS data acquisition was performed with an fNIRS brain imaging system (NIRScout 24-24, NIRx Medizintechnik GmbH, Berlin, Germany). This instrument has a laser near-infrared light source with four wavelengths at 685 nm, 780 nm, 808 nm, and 830 nm. The arrangements of the optode structure and fNIRS channels are shown in Figure 4 A, with 24 NIR sources, 24 NIR detectors, and 81 fNIRS channels. As can be seen in Figure 4 A, the fNIRS channels cover the whole brain cortex in order to analyze cortical connectivity, and it shows the experimental conditions used during fNIRS. fNIRS data were obtained during the pre-sequential finger-tapping task and post-sequential finger tapping task, which have the task block before and after tDCS according to the experimental design in Figure 4 B. The finger-tapping task paradigm starts with the first one-minute eye closed state to obtain a baseline fNIRS signal. The finger-tapping task was performed for 20 seconds with sequentially provided visual stimulation that permitted a finger press in the same position as the star mark ( ). And then, a subject took a rest for 20 seconds and repeated these blocks five times pre-and post-tDCS.

B. tDCS SYSTEM
As can be seen in Figure 4 A, the transcranial direct current stimulation (tDCS) electrodes were placed on the scalp, with the anode electrode overlying the right M1 region, which is the area related to the finger-tapping task of the lefthand, that was connected to the cathode electrode overlying the left M1 region. Stimulation was applied using a DC-STIMULATOR (NeuroConn, neuroCare Group, GmbH, Berlin, Germany). The DC current was initially increased in a ramp-like fashion over several seconds (∼10 s) until reaching 1 mA, and stimulation was maintained for a total of 30 min. DC currents were turned off slowly over a few seconds, out of the field of view of the patients, a procedure that does not elicit perceived sensations [34].

C. FUNCTIONAL CONNECTIVITY ANALYSIS
The sequential finger-tapping task was performed before and after tDCS as shown in Figure 4 B. The experimental block design for the finger tapping task consists of 20-sec task blocks alternating with 20-sec rest blocks to prevent finger fatigue. Five sets of alternating task and rest blocks were performed, and a total of 260 seconds of fNIRS data were acquired, including 1 minute of baseline fNIRS data for twenty-three participants and a total of 115 task blocks.  Estimation of the oxy-hemoglobin (HbO) signal from fNIRS was performed by using the NIRS-SPM toolbox [15] in MATLAB. Band-pass filter and normalization were used as the preprocessing methods for the raw HbO signals. The cutoff frequencies were set between 0.02 ∼ 0.1 Hz to remove long-term baseline drifts, high-frequency noise, and cardiac pulsations [35], [36]. And then, band-passed fNIRS data were normalized according to each task block. The fNIRS channels for each functional brain region were selected through the existing toolbox called fNIRS optodes' location decider (fOLD), which is a first-order approach to bring the advanced parcellation methods and meta-analyses acquired from functional magnetic resonance imaging to more precisely guide the selection of optode positions for fNIRS experiments [37]. The decided fNIRS channels were grouped as follows: MPF (CH 19, 77, 78 (CH 45, 46, 47), Rt. Pr (CH 51, 52, 53), and Occ (CH 51). In order to avoid signal distortion from differences in the number of fNIRS channels in each group, the fNIRS signals were processed with a narrow bandpass filter and normalization, then they were grouped through ensemble averaging.
Functional connectivity was estimated using the PLV [33], and group analysis for all task blocks of all participants was performed with the use of OptoNet II based on estimations using all trials and all epochs, as shown in Figure 5. The functional connectivity results obtained using the original version of OptoNet are shown in Figure 5 A, and were obtained based on fNIRS channels with manual thresholding (PLV > 0.8). As results for connectivity between functional brain regions cannot be analyzed for significance using the original OptoNet due to its fNIRS channel-based analysis and manual thresholding method, the results in Figure 5 A show all connectivity, even insignificant connectivity. Therefore, it is not easy to interpret these results, because there are so many connections among the fNIRS channels. Figure 5 B, C, and D show the results of connectivity analyses using the new OptoNet II. Results concisely presented connectivity between the functional brain regions and showed only connections that were significant using automatic thresholding with the random surrogate tests of 10,000 iterations. Figure 5 B represents the connectivity analysis performed before tDCS, which demonstrated several weak connections spread over wide brain areas between the frontal, temporal, M1, sensory, parietal, and occipital cortices. In particular, the most strong connection was seen between the occipital and parietal areas which can be interpreted as neural activities during the process of performing a motor task using visual cues [34]. Figure 5 C demonstrates the results of post-tDCS connectivity analysis that showed fewer significant connections. Instead, significant new connections were seen between M1, the sensory cortex, and parietal areas following tDCS. This is a result of the changes in the regions related to the motor task and the other associated regions because the brain networks around M1 and sensory cortices are increased after stimulation and the other networks not associated with motor control are decreased after stimulation. This is in line with previous studies, and it is thought that cortical networks begin to work efficiently and form essential connections during a motor task following tDCS [38]- [40]. However, to more clearly identify the difference before and after tDCS, the difference in connectivity before and after stimulation was estimated by using the new toolbox. Figure 5 D shows the differences in connectivity between Figure 5 C and B. As shown here, the significantly stronger connections following tDCS exist between the M1, temporal, parietal, and occipital regions. These are important connections related to the finger-tapping task with attention to a visual stimulus. The connections between the M1s of both hemispheres and adjacent cortical regions are considered to be the most important brain networks related to the motor task used in this experiment. By using OptoNet II, these networks were visualized by calculating the differences in connectivity strength between the pre-and post-tDCS conditions Therefore, the results of this analysis can support that changes in connectivity between functional brain regions related to the finger tapping task were enhanced following tDCS.

IV. DISCUSSION
In 2019, we released a MATLAB toolbox, named OptoNet, which can be used to analyze functional cortical networks for fNIRS. It was easy, and simple to use for plotting fNIRS signals and functional cortical connections according to fNIRS channels without any anatomical information. However, if there are many fNIRS channels, it might be difficult for OptoNet to find any given important connection due to so many connections.
In order to overcome these issues, we updated OptoNet functions. The first updated function is functional brain region-based analysis, which can easily save and load a brain region model after simply entering the number of fNIRS channels, and it can allow for analysis of the connectivity between functional brain region from fNIRS data. The second new capacity is analysis of connectivity differences, which is a function of loading two different saved connectivity results, then representing the different connectivity by task or state between them. We confirmed the clearer effect of tDCS through the connectivity difference between before and after tDCS stimulation. The last new function in OptoNet II is random surrogate tests-based auto-thresholding, which allows users to automatically identify only significant connections. The automatic threshold method shows objective meaningful connections by statistical significance through random comparisons of more than 10,000 iterations rather than just a cut off based in simple numerical comparison. Through these updates to OptoNet II, we are able to provide much better results reliability to OptoNet II users.
In this work, we developed a freely downloadable MATLAB toolbox, named OptoNet II (https://sites.google. com/site/dsucore/free/optonet) for the analysis of functional cortical connectivity. Any constructive comments and questions are always welcomed through our e-mail.
For future studies, we will update the OptoNet II toolbox to have even more user-friendly functions. More methods to evaluate functional connectivity networks to be added include directed transfer function (DTF), Granger causality, transfer entropy, partial directed coherence, and others as programming advances continue. We also plan to enable this toolbox not only for functional network analysis but also for brain activity analysis so that it can be used in all fNIRS research fields. YUN-HEE KIM is currently a Professor with the Department of Physical and Rehabilitation Medicine and serves as the Director for the Center for Prevention and Rehabilitation, Heart Vascular Stroke Institute, Samsung Medical Center, School of Medicine, Sungkyunkwan University, Seoul, South Korea. Her clinical specialty is neurological rehabilitation for patients with brain disease. She has authored more than 260 peer-reviewed articles and ten book chapters in the field of neurologic rehabilitation. Her current research interests include human neural plasticity using functional neuroimaging, noninvasive brain stimulation such as transcranial magnetic or direct current stimulations, and rehabilitation robotics. She is the President of the Korean Society for Neurorehabilitation and also holds other leadership positions in other Korean academic societies. VOLUME 9, 2021