On fractality of functional near-infrared spectroscopy signals: analysis and applications

Abstract. Significance: The human brain is a highly complex system with nonlinear, dynamic behavior. A majority of brain imaging studies employing functional near-infrared spectroscopy (fNIRS), however, have considered only the spatial domain and have ignored the temporal properties of fNIRS recordings. Methods capable of revealing nonlinearities in fNIRS recordings can provide new insights about how the brain functions. Aim: The temporal characteristics of fNIRS signals are explored by comprehensively investigating their fractal properties. Approach: Fractality of fNIRS signals is analyzed using scaled windowed variance (SWV), as well as using visibility graph (VG), a method which converts a given time series into a graph. Additionally, the fractality of fNIRS signals obtained under resting-state and task-based conditions is compared, and the application of fractality in differentiating brain states is demonstrated for the first time via various classification approaches. Results: Results from SWV analysis show the existence of high fractality in fNIRS recordings. It is shown that differences in the temporal characteristics of fNIRS signals related to task-based and resting-state conditions can be revealed via the VGs constructed for each case. Conclusions: fNIRS recordings, regardless of the experimental conditions, exhibit high fractality. Furthermore, VG-based metrics can be employed to differentiate rest and task-execution brain states.


Introduction
Functional near-infrared spectroscopy (fNIRS) is a promising noninvasive imaging technique that uses light in the near-infrared range to measure the local changes in the oxy-(Δ½HbO 2 ) and deoxygenated hemoglobin (Δ½HbR) concentrations associated with the underlying brain activities. 1 Compared to functional magnetic resonance imaging (fMRI), which is only sensitive to Δ½HbR, fNIRS provides additional information related to brain activity by also measuring Δ½HbO 2 . Additionally, fNIRS is portable and relatively less expensive compared to most neuroimaging modalities, enabling studying the brain function in more realistic settings. [2][3][4][5][6] Because of these advantages, fNIRS has been used in several cognitive and clinical neuroscience studies [7][8][9][10][11][12][13][14] as well as those involving brain-computer interfaces (BCIs). [15][16][17][18][19] A majority of fNIRS brain imaging studies have focused on the spatial domain and typically have ignored consideration of changes that occur in the temporal domain. Examples include *Address all correspondence to Laleh Najafizadeh, E-mail: laleh.najafizadeh@rutgers.edu localization studies, [20][21][22] where the aim is to identify brain's activation patterns in response to specific stimuli, and connectivity studies (functional or effective), where the focus is on investigating the functional interactions among brain regions, either when the brain is at rest or is engaged in performing a particular task. [23][24][25][26][27] However, it is now well known that the brain is highly dynamic, [28][29][30][31][32] and therefore, to gain a more comprehensive picture about its function, methods capable of extracting temporal information in brain recordings are required.
As compared to the spatial domain, a much smaller number of fNIRS studies exists that have considered the temporal domain for the analysis. [33][34][35][36][37][38][39][40] For example, in Ref. 33, by applying the Higuchi fractal dimension algorithm, 41 it is shown that fNIRS signals have high degree of complexity. The wavelet transform is applied to fNIRS signals, and it is shown that the wavelet coefficients can be used to train a classifier. In Refs. 38-40, entropy has been used to assess the complexity of fNIRS signals in patient groups (such as those with Alzheimer's disease, attention-deficit hyperactivity disorder, and traumatic brain injury), demonstrating that it carries information that can be related to diseases. All these studies suggest that there exists information relevant to the underlying brain activity in the complex characteristics of fNIRS signals.
In this paper, using visibility graph (VG), we present an approach for revealing the fractal properties of fNIRS time series. VG is a recently introduced method, which maps a time series into a graph (called a VG). As will be discussed, the topological properties of the constructed graph are related to the fractality and complexity of the time series. 42,43 Compared to conventional fractal analysis approaches, 42 VG is computationally less complex and has been used in various studies. [44][45][46][47][48][49] For example, using electrocardiogram, Jiang et al. showed that employing VG analysis can reveal the dynamical changes caused by mediation training, manifested as regular heartbeat, which is closely related to the adjustment of the autonomous neural system. 44 Zhu et al. applied a VG-based approach for alcoholism identification, showing that this approach is promising in separating alcoholic subjects from controlled drinkers. 48 In Ref. 47, it was shown that VG applied to electroencephalography (EEG) signals can provide features that can distinguish children with autism from nonautistic children. In Ref. 49, we have shown that the temporal characteristics of calcium recordings in GCaMP6 mice extracted by VG carries discriminatory information that can be utilized to decode behavior.
It is important to note here the difference between VG and the graph theoretical-based methods used commonly in functional connectivity studies. 50,51 In typical functional connectivity studies, graphs are constructed in the spatial domain, i.e., nodes in the graph correspond to the location of channels or voxels, and links between two nodes are formed based on the statistical similarity of the time series associated with the two nodes, quantified by measures, such as correlation. On the other hand, as will be discussed in Sec. 2, in VG, the nodes correspond to the time points in the time series, and the links are formed based on natural visibility between the time points ( Fig. 1). Once the graph is formed for each time series, graph measures can be extracted to represent different properties of the time series.
In this paper, we use VG to study the fractality of fNIRS time series under two conditions: when the brain is at rest and when it is engaged in the task. fNIRS time series are recorded from nine healthy male subjects at two resting-state conditions and two task conditions. VGs are constructed from each time series for each channel and each condition. The power of scale-freeness of visibility graph (PSVG) is then extracted and compared across conditions.
To the best of our knowledge, this is the first study to employ VG to reveal temporal characteristics of fNIRS-recorded time series, demonstrating its feasibility in identifying features in fNIRS recordings that can be used to gain new insights about the brain function.
The rest of this paper is organized as follows. Section 2 describes the methods used for analysis in this study. The details of the experimental setup are given in Sec. 3. Section 4 presents the results, and finally, some discussions are provided in Sec. 5.

Methods
Fractal analysis of time series provides an interesting opportunity to study their temporal structure in terms of self-similarity, power-law scaling relationship, and scale-invariance. 52 Since its introduction, fractal analysis has been used in several fields, such as physics, geography, biology, and psychology 53 to reveal properties in time series that will not be visible through conventional analysis. Recently, fractal analysis has also been applied to recordings obtained through functional brain imaging studies 42,[54][55][56] to better characterize the observed temporal fluctuations of the signals or to differentiate diseased group from healthy population with an improved accuracy. For instance, using fractal analysis of fMRI data, Sokunbi et al. showed that the hemodynamic response measured from patients with Schizophrenia presented larger complexity, compared to that measured from healthy controls. 57 He showed that the fractal properties of the fMRI signal altered with the changes of brain functional state. 58 Using fNIRS, Khoa and Nakagawa showed that the complex characteristics of the signals recorded during physical motion and imaginary motion of the right arm were different, and these differences can be potentially used in BCIs. 33 Lei et al. showed that the power spectra of both EEG and fMRI signals follow the power-law distribution, and scale-free brain activity can be characterized by robust temporal structures. 59 There exist several methods to estimate the fractality of a time series, 52 and the choice of the proper method has to be made based on the properties of the time series, and whether it is stationary or nonstationary. In this paper, the fractal properties of fNIRS time series will be evaluated using two approaches: scaled windowed variance (SWV) analysis and VG. These methods are described here.
In the following, we denote an N-point time series (e.g., recorded by a given fNIRS channel) with x ¼ fx i g N i¼1 .

Scaled Windowed Variance Analysis
We first evaluate the fractal dimension of fNIRS time series using conventional methods. As it is well known that physiological signals and brain activities are nonstationary, 52,60 SWV analysis method 52,61 is used here to estimate the fractal dimension (defined below) of fNIRS-recorded time series.
To estimate the fractal dimension of a time series using SWV, the time series is partitioned into nonoverlapping segments of size n. Ifμ represents the mean of the segment, the standard deviation for each segment,σ n , is computed as follows: 52 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 1 4 0σ This measure is computed for all segments and then is averaged to obtain σ n . The procedure is repeated for different window sizes. For nonstationary fractal time series, the windowed mean standard deviation, σ n , and its corresponding window size n follow a power law relation given by 52 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 1 1 where ¼ d represents equal in distribution, H is the Hurst coefficient that can be obtained by calculating the slope of the least-squares linear regression line of logðσ n Þ versus logðnÞ, and p is a proper prefactor. 52 The fractal dimension, D, is linearly related to H as D ¼ 2 − H. 62

Visibility Graph
VG is a recently introduced method for studying the fractal properties of time series. 43 The VG associated with a given time series x of N points is constructed as follows (see Fig. 1). Each time point in x is considered as a node in the graph (i.e., for an N-point time series, the graph will have N nodes). The link between two nodes is formed if the nodes are considered to be "naturally visible." That is, in the graph, there will be a link between nodes h and l (corresponding to time points t h and t l in the time series), if and only if, for any node p (t h < t p < t l ), the following condition holds: 43 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 1 5 The VG is expressed using the adjacency matrix, an N × N symmetric matrix A N×N ¼ ½a hl , where a hl ¼ 1 when nodes h and l are connected and a hl ¼ 0, otherwise. 42 The graph that is constructed through VG reveals the dynamic properties of the time series in unique ways. For example, periodic signals result in regular graphs, and fractal time series result in scale-free networks. 42,43,45,63 Scale-free corresponds to the property of the graph that, independent of the number of nodes, its degree distribution PðkÞ has a power-law tail, where the tail exponent obeys the power law, i.e., E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 8 6 PðkÞ ∼ k −γ : In Eq. (4), k represents the degree of the node (i.e., the number of links connected to a node), PðkÞ denotes the degree distribution (i.e., the fraction of nodes with degree k), and γ denotes the PSVG. It has been proven that the PSVG is indicative of the fractality of the time series. 64

Experimental Procedure
In this section, we describe the experimental setup and the preprocessing steps used to remove artifacts from fNIRS signals. Figure 2 shows the experimental setup and the location of optodes.  65 Participants were asked to click only when the Go symbol is presented. Recordings for the two resting-state and the two tasks sessions were all completed in one setting. Changes in the optical signals were recorded at two wavelengths (685 and 830 nm) at the sampling rate of 12.5 Hz, using NIRScout system (NIRx Medical Technologies, LLC). Note that although the maximum sampling rate of NIRScout system is 62.5 Hz, in this experiment, 12.5 Hz was used as the practical sampling rate due to the time-division multiplexing on the illuminated light sources to avoid cross talk. A total of 14 channels were used (8 light sources and 8 detectors arrangement). For each channel, the distance between the light source and the light detector was 3 cm. No short distance probes 66 were used. Optodes were placed to cover the prefrontal and motor cortices, where brain activities are expected for these tasks according to prior studies. 65,67,68

Data Preprocessing
In this section, we describe the artifact rejection, bandpass filtering, and hemodynamic signal conversion procedures.
Motion-related artifacts inevitably exist in fNIRS-recorded time series, and if not removed, can negatively impact the outcome of fractal analysis. Therefore, artifact rejection procedure was first performed. 69 Two types of artifacts, "spikes" and "discontinuities" (jumps), were particularly considered, as they could impact the construction of VG. For each channel, short intervals containing spikes were identified by visual inspection and then replaced by the average of the signals of the same interval from two neighboring channels. Discontinuities were detected by examining the difference in signal values at every two successive time points. For every time point pair, discontinuity was detected if this difference becomes larger than four times of the standard deviation of the time series. In such cases, the discontinuity was eliminated by subtracting the difference from the values for the point after the jump. 70 Next, optical signals were band-pass filtered using a fourth-order band-pass Butterworth filter in the range of 0.01 to 0.2 Hz to exclude the low-frequency drift and the interference caused by physiological sources, such as heart beat and respiration. 71,72 The filtered signals were converted into Δ½HbO 2 and Δ½HbR according to the modified Beer-Lambert's law, based on the following equations: 73 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 2 1 2 ln I λ 1 Note that Δ½HbO 2 and Δ½HbR are the unknown parameters here and will be found by solving this system of two linear equations. In Eq. (5), I λ i and I baseline;λ i (i ¼ 1;2) are the optical intensities measured at the detector location at wavelength λ i during experimental block and during preblock baseline period, respectively, x is the distance between the light source and the light detector, DPF λ i is the differential pathlength factor, and ϵ HbO 2 ;λ i and ϵ HbR;λ i are the extinction coefficients of HbO 2 and HbR at wavelength λ i , respectively. Here, I baseline;λ i for λ i is obtained as the averaged recorded light intensity, associated with wavelength λ i , during the 10-min EC resting-state session. The artifact rejection, band-pass filtering, and data conversion procedures were performed using the nirsLAB toolbox. 70,74 For the RT and GNG tasks, time series were segmented according to the border of the blocks (170 s for each block). The same duration was used to segment the signals from the resting-state recordings into three nonoverlapping blocks. This procedure for all subjects resulted in a total of 108 fourteen-channel time series, each with 2125 data points.

Results
In this section, we present the results of fractal analysis from the SWV and VG approaches.

Results from Scaled Windowed Variance Analysis
Fractality of fNIRS recordings, for each channel, subject, condition, and block, was first examined using SWV analysis, which is an appropriate method for nonstationary fractal time series. 52 Figure 3(a) shows the obtained mean standard deviation [σ x ðnÞ] as a function of the window size (n) in a double logarithmic plot for a given time series. As can be seen, for the window size ranging from 2 2 to 2 7 , the data follow a linear trend. By definition, the slope of this trend equals the Hurst exponent, H, from which the fractal dimension D can be estimated from D ¼ 2 − H. 52 As a measure of self-similarity, the value of D ranges between 1.0 and 2.0. D for a random time series is 1.5. For a time series with high correlation or memory over time, the extracted D value is near 1. 62 The distribution of estimated D values for all preprocessed time series (resting-state and task) is shown in Fig. 3(b). As can be seen, the estimated D value for all time series is less than 1.3 with a mean of 1.15. This result confirms the existence of high degree of fractality (positive correlation) in fNIRS time series, and hence, VG analysis can be applied to these signals.

Results from Visibility Graph Analysis
We first examine that the temporal structure of the fNIRS signals can be characterized using the PSVG measure. To achieve this goal, the VG was constructed for two cases: for a representative preprocessed fNIRS time series and for its randomly shuffled version (i.e., the order of appearance of data points in time was randomly shuffled). The PSVG was then estimated from VG for each case. Figure 4(a) shows an example of a time series (a Δ½HbO 2 signal associated with one EC block) and its associated power-law tail of the degree distribution. The shuffled version of this time series and its associated power-law tail of the degree distribution are shown in Fig. 4(b). It can be seen that the degree distribution in the two cases decays at different rates. The estimated PSVG for the original time series [ Fig. 4(a)], PSVG orig , is 3.512. For the shuffled version, random shuffling was performed for 100 rounds, and for each round, the power of scale-freeness, PSVG shuf , was estimated. The averaged PSVG shuf is obtained as 3.996 [standard deviation is 0.120, range [3.787 to 4.283], refer to Fig. 4(c)]. It is observed that although both the original time series and its randomly shuffled version have the same distribution in terms of amplitude of data points, they have different PSVG values. This result implies that the temporal structure of the time series, rather than the distribution of amplitude of data points, can be characterized from the PSVG values.
Next, the VGs were then constructed for the time series associated with each channel, subject, condition, and block separately. As an example, Figs. 5(a) and 5(b) show two representative Δ½HbO 2 signals recorded from one subject under EC and GNG blocks, respectively. The waveforms with different colors represent the Δ½HbO 2 time series recorded from different channels. Figures 5(c) and 5(d) illustrate two adjacency matrices of VGs associated with EC and GNG from channel 1, respectively. It can be seen that there exist differences in the number of links of the two matrices, indicative of differences in temporal structure of the two time series.
To further quantify these differences, we extracted the graph measure of PSVG for each case. Figure 6(a) presents the degree distribution patterns averaged across subjects for each of the 14 channels. The colors correspond to different conditions (EC, EO, RT, and GNG). The zoomed-in power-law tail of the distributions is shown in Fig. 6(b). As can be seen, for majority of channels, the slopes are different for resting-state and task-based time-series.
A least-square regression line was fitted to the power-law tail and from it, the PSVG (equal to the negative of the slope) was calculated (γ ¼ − log½PðkÞ logðkÞ ). The mean and standard deviation of the estimated PSVG values associated with each condition and each channel are presented as bar plots in Fig. 6(c).
To evaluate if the observed differences between PSVG values across conditions are statistically significant, we first performed a one-way repeated analysis of variance (ANOVA). Results are summarized in Table 1. The results from ANOVA with majority of channels showing p < 0.05 further motivated us to perform paired-sample student's t-tests for each channel across each pair of conditions (EO-EC, EO-RT, EO-GNG, EC-RT, EC-GNG, RT-GNG). For each channel, the pairs of conditions that result in significant difference were identified and labeled using "star" notations in Fig. 6(c). The results of this statistical test are also summarized in  Table 2. It is shown clearly that for all subjects, the PSVG values associated with resting states (EC and EO) are larger than those associated with task-based (RT and GNG) conditions. Using this result, we hypothesize that the PSVG values of VGs can be used to distinguish brain states. In the following, we perform classification to test this hypothesis.

Classification Results
Binary classification was performed to evaluate the feasibility of using PSVG in distinguishing brain states. The individual multichannel PSVG values associated with EC and EO conditions were pooled to form the resting-state sample data, and those PSVG values associated with RT and GNG conditions were pooled to form the task-execution sample data. Before pooling, the individual PSVG values were normalized across conditions, so their l 2 -norm was equal to 1 for each channel. This procedure left 18 data for each binary condition. For the classifiers, we chose kNN, linear and nonlinear support vector machine (SVM), linear discriminant analysis (LDA), and quadratic discriminant analysis (QDA). We note that previous studies have used artificial neural networks (ANN) for classification of fNIRS data, 75 but considering the small sample size in this study, we did not implement ANN.
Given the small sample size here, leave-two-out-cross-validation procedure was used, in which the classification was repeated   The results for accuracy, specificity, and sensitivity are summarized in Table 3. The best result for kNN is achieved for k ¼ 5. The linear SVM performs better than nonlinear SVM. For DA, the performance of LDA is better than QDA. It is observed that the nonlinear classifiers (QDA and nonlinear SVM) did not improve the accuracy of classification due to the relatively small sample size. Overall, it can be seen that kNN, SVM, and LDA have shown significantly better accuracy than that of the chance level (50%). These results reveal that the PSVGs of hemodynamic signals measured from both prefrontal and motor cortices carry discriminatory information for differentiating resting-state and task-execution conditions.

Discussions and Conclusions
The human brain is a highly complex system with nonlinear dynamical behavior. 76 As such, methods such as fractal analysis, which can reveal nonlinearities in the recordings associated with brain activities, can serve as an important complementary approach to the commonly used  methods (e.g., used for functional localization or functional connectivity) in order to gain a more comprehensive knowledge about how the brain functions.
In this study, we investigated the fractality of fNIRS-recorded time series obtained during task execution or when the brain is at rest using VG. Fractal dimension for all conditions and recordings was first estimated through SWV analysis. It was shown that fNIRS recordings, regardless of the experimental conditions, exhibit high fractality. This result is in line with previous studies that used EEG or fMRI to monitor brain function. 42,77,78 Next, VGs were constructed for all fNIRS time series recordings and their corresponding PSVG values were calculated. Results showed that for most channels, the difference in PSVG values for cases when the brain is at rest and when the brain is engaged in executing tasks is statistically significant. To the best of our knowledge, this is the first study that explores the possibility of employing VG in differentiating rest and task-execution brain states.
The capability of VG-based metrics in differentiating brain states was further examined by performing classification using the PSVG values estimated during rest or task execution as features to the classifier. The feature dimension was equal to the number of channels. A wide range of classifiers was used and although the number of training and testing samples in this study was small, a reasonably good accuracy was obtained.
One application of VG-based metrics being used as features would be to use them as biomarkers for diagnosing brain-related disorders. [38][39][40] Due to the instrumental advantages of fNIRS technology (e.g., portability and low cost), fNIRS has great potential for clinical settings. As such, finding features in fNIRS recordings that can serve as biomarkers for diagnosis are of great interest to clinicians. The results of this work show the capability of VG-based metrics in finding features that are based on nonlinear properties of fNIRS recordings, which could have applications in clinical settings, for example, in cases, where disease-specific information exist in the nonlinear temporal characteristics of the recordings.
In future studies, we plan to investigate how other metrics of VG, such as closeness and betweenness centrality, clustering coefficient, and transitivity, 49,63 can characterize the temporal structures of fNIRS recordings. Additionally, we plan to extend VG to multilayer VG, where relations among time series of different channels can also be quantitatively characterized, 79 providing more interesting information about how the brain functions.

Disclosures
The authors declare that there are no conflicts of interest related to this article.