*Article* **Sample Entropy Combined with the K-Means Clustering Algorithm Reveals Six Functional Networks of the Brain**

#### **Yanbing Jia <sup>1</sup> and Huaguang Gu 2,\***


Received: 30 September 2019; Accepted: 22 November 2019; Published: 26 November 2019 -

**Abstract:** Identifying brain regions contained in brain functional networks and functions of brain functional networks is of great significance in understanding the complexity of the human brain. The 160 regions of interest (ROIs) in the human brain determined by the Dosenbach's template have been divided into six functional networks with different functions. In the present paper, the complexity of the human brain is characterized by the sample entropy (SampEn) of dynamic functional connectivity (FC) which is obtained by analyzing the resting-state functional magnetic resonance imaging (fMRI) data acquired from healthy participants. The 160 ROIs are clustered into six clusters by applying the K-means clustering algorithm to the SampEn of dynamic FC as well as the static FC which is also obtained by analyzing the resting-state fMRI data. The six clusters obtained from the SampEn of dynamic FC and the static FC show very high overlap and consistency ratios with the six functional networks. Furthermore, for four of six clusters, the overlap ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC, and for five of six clusters, the consistency ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC. The results show that the combination of machine learning methods and the FC obtained using the blood oxygenation level-dependent (BOLD) signals can identify the functional networks of the human brain, and nonlinear dynamic characteristics of the FC are more effective than the static characteristics of the FC in identifying brain functional networks and the complexity of the human brain.

**Keywords:** sample entropy; brain functional networks; complexity; dynamic functional connectivity; static functional connectivity; K-means clustering algorithm

#### **1. Introduction**

The human brain shows complex spatiotemporal behaviors when executing physiological functions. Characterizing dynamics of the complex spatiotemporal behaviors is of great significance in understanding the human brain. Since blood oxygenation level-dependent (BOLD) signals of different brain regions can be measured by the functional magnetic resonance imaging (fMRI) technique at high spatial and temporal resolutions, BOLD signals have been widely used to characterize dynamics of the spatiotemporal behaviors of the human brain [1,2]. For instance, the temporal correlation in BOLD signals of two distinct brain regions is commonly employed to describe the functional connectivity (FC) between them [3]. A positive and strong temporal correlation corresponds to a strong FC, and some brain regions with strong FCs among them constitute a brain functional network [4–6]. Alterations of some FCs in a brain functional network are often associated with brain disorder, such

as schizophrenia [7], major depression [8], autism [9], Alzheimer's Disease [10], and attention deficit hyperactivity disorder [11]. For example, Cheng et al. evaluated the FC between different brain regions in subjects with autism and found a key system in the middle temporal gyrus with reduced FC and a key system in the precuneus with reduced FC [12].

In most previous research on FC, only one correlation coefficient is acquired using entire BOLD signals of two distinct brain regions. The one correlation coefficient is called the static FC between the two brain regions. Recently, to understand dynamics of the spatiotemporal behaviors of the human brain more deeply, some researchers acquired a sequence of correlation coefficients by applying the sliding-window approach to BOLD signals of two distinct brain regions [13–23]. These correlation coefficients form a time series which is called the dynamic FC between the two brain regions. The dynamic FC exhibits complex characteristics which are effective in describing properties of the brain functional networks of patients with brain disorder. For instance, in one of our recent studies, complex characteristics of dynamic FC were described by sample entropy (SampEn), and the effects of schizophrenia on such complex characteristics were investigated. It was shown that the visual cortex of the patients with schizophrenia exhibited significantly higher SampEn than that of the healthy controls [24]. As introduced above, both the static FC and the SampEn of dynamic FC are effective in describing properties of the brain functional networks of patients with brain disorder. However, the effectivenesses of the static FC and the dynamic FC have not been compared directly.

Studies on the static FC or the dynamic FC are often carried out by first extracting BOLD signals of different brain regions and then evaluating the static or the dynamic FC between different brain regions for further analysis. Different brain regions are often determined by a brain template, such as the Dosenbach's template [25]. The Dosenbach's template includes 160 regions of interest (ROIs) determined by a sequence of meta-analyses of task-based fMRI studies which cover much of the human brain [25]. Furthermore, the 160 ROIs can be separated into six functional networks including the default, the frontal-parietal, the cingulo-opercular, the sensorimotor, the occipital, and the cerebellum networks, which were identified by performing modularity optimization on the average FC matrix across a large cohort of healthy subjects [25]. The six functional networks have been used in predicting brain maturity across development [25,26], parcellating cortical or subcortical regions [27], examining the influence of temporal properties of BOLD signals on FC [28] and so on. For instance, Zhong et al. parcellated the hippocampus based on the FC, and showed that both the left and right hippocampus were divided into three subregions exhibiting different FC profiles with the six functional networks [27]. However, machine learning algorithms have not been used to identify the six functional networks.

The K-means clustering algorithm is one of the unsupervised learning algorithms [29]. Since the K-means clustering algorithm can cluster different observations into different clusters in a simple and easy way, it has been widely used in fMRI studies [30–38]. For instance, Fan et al. used the K-means clustering algorithm to parcellate the thalamus based on the static FC and found that the thalamus could be divided into seven symmetric thalamic clusters [36]. Park et al. parcellated the primary and secondary visual cortices (V1 and V2) into several subregions by applying the K-means clustering algorithm to the static FC and found that V1 and V2 could be separated into anterior and posterior subregions [38].

The present study intends to cluster the Dosenbach's 160 ROIs into six clusters by applying the K-means clustering algorithm to the static FC and the SampEn of dynamic FC, to analyze the overlap and consistency between the six clusters and the six functional networks, and to compare the effectivenesses of the static FC and the dynamic FC. It is shown that applying the K-means clustering algorithm to FC is feasible to identify the six functional networks, and the SampEn of dynamic FC is more effective than the static FC as the six clusters obtained from the SampEn of dynamic FC show higher overlap and consistency ratios with the six functional networks.

This paper is organized as follows. The experiments and methods are presented in Section 2. The cluster results for the static FC and the SampEn of dynamic FC and the comparisons between them are shown in Section 3. The conclusion and discussion are described in Section 4. Some supplementary tables are presented in the appendix.

#### **2. Experiments and Methods**

#### *2.1. Participants*

FMRI data for this study were acquired at Olin Neuropsychiatry Research Center and have been made publicly available http://fcon\_1000.projects.nitrc.org/indi/abide. The data were acquired from 31 healthy participants (18 males and 13 females) over the age range 18–30 years. This sample was retained after applying criteria for head motion, from a total of 35 healthy participants. Informed consent was obtained from all participants in accordance with Olin Neuropsychiatry Research Center Institutional Review Board oversight.

#### *2.2. Data Acquisition and Preprocessing*

BOLD signals are extracted from three-dimensional functional images collected on a Siemens 3T MRI scanner with the following parameters: repetition time (TR), 475 ms; echo time, 30 ms; field of view, 240 × <sup>240</sup> mm2; slices, 48; slice thickness, 3 mm; flip angle, 60◦. During the data collection, all participants were instructed to rest but not fall asleep. For each participant, 947 three-dimensional functional images were collected.

The functional images are preprocessed using SPM8 and DPABI softwares [39,40]. Firstly, the first 4 images are discarded to reduce the negative effects of scanner's stabilization on the analysis results. Secondly, the images are corrected for time delay in slice acquisition and rigid-body head motion. Thirdly, several confounding factors are regressed out from the images, including 6 head motion parameters and the cerebrospinal, the white matter, and the global brain signals. Fourthly, temporal band-pass filtering (0.01–0.08 Hz) of the images are performed to reduce the negative effects of low-frequency drift and high-frequency physiological noise on the analysis results. Fifthly, the images are spatially normalized to the Montreal Neurological Institute space and are resampled to voxels of size 3 × <sup>3</sup> × <sup>3</sup> mm3. Sixthly, the images are smoothed with a Gaussian kernel of 8 mm full-width at half-maximum. Finally, the BOLD signal of each voxel is extracted from the functional images.

#### *2.3. The Dosenbach's Template and the 6 Functional Networks*

One hundred and sixty regions of interest (ROIs) are selected based on the Dosenbach's template [25]. The centroid of each ROI is derived from a sequence of meta-analyses of task-based fMRI studies (Figure 1a). The radius of each ROI equals 5 mm (Figure 1a). The name and the sequential number of each ROI can be found in Table A1 in Appendix A. The 160 ROIs can further be grouped into 6 functional networks, including the default, the frontal-parietal, the cingulo-opercular, the sensorimotor, the occipital, and the cerebellum networks (Figure 1a). The name and the sequential number of each ROI in each functional network can be found in the first and second columns of Tables A2–A7 in Appendix A.

Based on the 6 functional networks, an adjacent matrix can be generated [36,41,42]. The adjacent matrix is labeled as

$$A = \begin{bmatrix} a\_{1,1} & \cdots & a\_{1,160} \\ \vdots & \ddots & \vdots \\ a\_{160,1} & \cdots & a\_{160,160} \end{bmatrix} \tag{1}$$

Each of the elements on the main diagonal of *A* is 1. Other elements of *A* are defined as follows: *ai*,*<sup>j</sup>* = 1 if the *i*th ROI and the *j*th ROI are contained in the same functional network and *ai*,*<sup>j</sup>* = 0 otherwise (*i*, *j* = 1, 2, . . . , 160) (Figure 1b).

**Figure 1.** (**a**) One hundred and sixty regions of interest (ROIs) are shown on a surface rendering of the brain. ROIs in different functional networks are shown in different colors. (**b**) The adjacent matrix *A* of 160 ROIs in 6 functional networks.

#### *2.4. The Static FC and the Dynamic FC*

The BOLD signal of each ROI is extracted by averaging the BOLD signals over all voxels in this ROI. Then both the static FC and the dynamic FC are evaluated (Figure 2).

**Figure 2.** The static functional connectivity (FC) matrix *B* and the SampEn matrix *E* obtained from the BOLD signals of 160 ROIs. The matrices *B* and *E* are used to cluster the 160 ROIs into 6 clusters by the K-means clustering algorithm.

The static FC between each pair of ROIs is assessed by a Pearson correlation coefficient. For each of the 31 participants, after the static FC between each pair of ROIs is evaluated, a static FC matrix of size 160 × 160 is obtained (Figure 2), which is labeled as

$$B = \begin{bmatrix} b\_{1,1} & \cdots & b\_{1,160} \\ \vdots & \ddots & \vdots \\ b\_{160,1} & \cdots & b\_{160,160} \end{bmatrix} = \begin{bmatrix} B\_1 \\ \vdots \\ B\_{160} \end{bmatrix} . \tag{2}$$

The *i*th row *Bi* represents the static FC between the *i*th ROI and all the other ROIs (*i* = 1, 2, ... , 160). The matrix *B* is used to cluster the 160 ROIs into 6 clusters.

Dynamic FC is assessed by the sliding-window approach. Specifically, a tapered window is created by convolving a rectangle window (size = 20 TRs = 9.5 s) with a Gaussian curve (standard deviation = 3 TRs) [14,15,23]. The window is used to extract BOLD signals in a step of 1 TR, leading to 923 time windows per subject (Figure 2). For the *k*th time window (*k* = 1, 2, ... , 923), a Pearson correlation coefficient is used to evaluate the FC between each pair of ROIs and thus a FC matrix of size 160 × 160, which is labeled as

$$D\_k = \begin{bmatrix} d\_{1,1,k} & \cdots & d\_{1,160,k} \\ \vdots & \ddots & \vdots \\ d\_{160,1,k} & \cdots & d\_{160,160,k} \end{bmatrix} \tag{3}$$

which is obtained for each subject (Figure 2). As *k* increases from 1 to 923, *di*,*j*,*<sup>k</sup>* forms a time series (*i*, *j* = 1, 2, ... , 160), which represents the temporal evolution of the FC between the *i*th and *j*th ROIs and is named as the dynamic FC (Figure 2). Since previous studies showed that the window of size 20 TRs captures more transient patterns in dynamic FC [23], the window size is fixed at 20 TRs throughout the study.

#### *2.5. SampEn of a Dynamic FC Time Series*

For each dynamic FC time series, *di*,*j*(*i*, *j* = 1, 2, ... , 160, *i* = *j*), the SampEn is calculated. For convenience, time series *di*,*<sup>j</sup>* is denoted by **x** = (*x*1, *x*2, ... , *xN*)(*N* = 923). SampEn of **x** is computed as follows [24,43–46].

Firstly, constructing embedding vectors **vi** = (*xi*, *xi*+1, ... , *xi*+*m*−1), in which *m* stands for the dimension of **vi**(1 ≤ *i* ≤ *N* − *m* + 1).

Secondly, define

$$\mathbf{C}\_{i}^{m} = \frac{1}{N - m} \sum\_{j=1, j \neq i}^{N-m+1} \Theta(r - \|\mathbf{v\_{i}} - \mathbf{v\_{j}}\|). \tag{4}$$

*r* stands for a tolerance value which is defined as *r* = *ε* · *σ***x**, where *ε* is a small parameter and *σ***<sup>x</sup>** is the standard deviation of **x**. Θ(·), the Heaviside function, which is defined as

$$\Theta(\mathbf{x}) = \begin{cases} \ 0, & \mathbf{x} < \mathbf{0}; \\\ 1, & \mathbf{x} \ge \mathbf{0}. \end{cases} \tag{5}$$

· represents the Chebyshev distance, i.e.,

$$\|\mathbf{v\_i} - \mathbf{v\_j}\| = \max(|\mathbf{x\_i} - \mathbf{x\_j}|, |\mathbf{x\_{i+1}} - \mathbf{x\_{j+1}}|, \dots, |\mathbf{x\_{i+m-1}} - \mathbf{x\_{j+m-1}}|).\tag{6}$$

Similarly, define

$$C\_{i}^{m+1} = \frac{1}{N - m - 1} \sum\_{j=1, j \neq i}^{N-m} \Theta(r - ||\mathbf{v\_i} - \mathbf{v\_j}||). \tag{7}$$

*Entropy* **2019**, *21*, 1156

Thirdly, in view of Equations (4) and (7), we define

$$\mathcal{U}^{\text{pr}} = \frac{1}{N - m + 1} \sum\_{i=1}^{N-m+1} \mathbb{C}\_i^{\text{pr}} \, , \tag{8}$$

and

$$
\mathcal{U}^{m+1} = \frac{1}{N-m} \sum\_{i=1}^{N-m} \mathbb{C}\_i^{m+1}.\tag{9}
$$

Finally, calculate SampEn of **x** as

$$\text{SampEn} = -\ln \frac{\mathcal{U}^{m+1}}{\mathcal{U}^m}.\tag{10}$$

The value of SampEn is not less than 0, and a larger value of SampEn means more complexity [47]. Similar to our previous study [24,43], *m* and *ε* are fixed at 2 and 0.2, respectively.

In addition, because *di*,*i*,*<sup>k</sup>* = 1(*i* = 1, 2, ... , 160, *k* = 1, 2, ... , 923), the SampEn of *di*,*<sup>i</sup>* equals 0 (*i* = 1, 2, ... , 160). Thus, for each participant, a SampEn matrix of size 160 × 160 is obtained (Figure 2). The SampEn matrix is labeled as

$$E = \begin{bmatrix} \varepsilon\_{1,1} & \cdots & \varepsilon\_{1,160} \\ \vdots & \ddots & \vdots \\ \varepsilon\_{160,1} & \cdots & \varepsilon\_{160,160} \end{bmatrix} = \begin{bmatrix} E\_1 \\ \vdots \\ E\_{160} \end{bmatrix}. \tag{11}$$

The element *ei*,*<sup>j</sup>* represents the SampEn of dynamic FC between the *i*th ROI and *j*th ROI (*i*, *j* = 1, 2, ... , 160). *ei*,*<sup>i</sup>* equals 0 (*i* = 1, 2, ... , 160). The matrix *E* is used to cluster the 160 ROIs into 6 clusters.

*2.6. Clustering ROIs into 6 Clusters by Applying the K-Means Clustering Algorithm to the Static FC Matrix*

For each of the 31 participants, there exists a static FC matrix *B* of size 160 × 160. The *i*th (1 ≤ *i* ≤ 160) row *Bi* = (*bi*,1, *bi*,2,..., *bi*,160) represents the static FC between the *i*th ROI and all the other ROIs.

The K-means clustering algorithm is commonly used to cluster different observations into different clusters based on the distance between these observations [29]. In the present paper, the K-means clustering algorithm is applied to the matrix *B* to cluster 160 ROIs into 6 clusters. Procedures of the algorithm are briefly described as follows.

First, select 6 rows from the matrix *B* and use these 6 rows as initial cluster centroids.

Secondly, calculate the squared Euclidean distance between each row and each initial cluster centroid, and then assign each row to the cluster with the closest centroid.

Thirdly, when all rows have been assigned, calculate the average of the rows in each cluster to obtain 6 new cluster centroids.

Finally, repeat the second and the third steps until the centroids no longer change.

The algorithm generates 6 clusters, and each cluster is composed of different rows of the matrix *B* (or of different ROIs). Based on the 6 clusters, an individual adjacent matrix of size 160 × 160 is generated [36,41,42]. The individual adjacent matrix is labeled as

$$F = \begin{bmatrix} f\_{1,1} & \cdots & f\_{1,160} \\ \vdots & \ddots & \vdots \\ f\_{160,1} & \cdots & f\_{160,160} \end{bmatrix} \tag{12}$$

Each of the elements on the main diagonal of *F* is 1, and other elements of *F* are defined as follows: *fi*,*<sup>j</sup>* = 1 if the *i*th ROI and the *j*th ROI are contained in the same cluster and *fi*,*<sup>j</sup>* = 0 otherwise.

Since the study includes 31 participants, 31 individual adjacent matrices are obtained. A group adjacent matrix of size 160 × 160 is obtained by averaging 31 individual adjacent matrices. The group adjacent matrix is labeled as

$$\mathbf{G} = \begin{bmatrix} \mathcal{S}1, 1 & \cdots & \mathcal{S}1, 160 \\ \vdots & \ddots & \vdots \\ \mathcal{S}160, 1 & \cdots & \mathcal{S}160, 160 \end{bmatrix} \cdot \tag{13}$$

The K-means clustering algorithm is further applied to the matrix *G* to obtain the group cluster result [36,41,42] and the 6 clusters of the group cluster result are compared with the 6 functional networks shown in Figure 1a.

The detailed clustering procedure is performed by MATLAB software (MATLAB R2014b). Considering that the K-means clustering algorithm is sensitive to the initial cluster centroids, we repeat each clustering procedure 500 times, and the cluster result with the lowest within-cluster distance is adopted.

#### *2.7. Clustering ROIs into 6 Clusters by Applying the K-Means Clustering Algorithm to the SampEn Matrix*

The procedures described in Section 2.6 are also applied to the SampEn matrix *E*, and 6 clusters are obtained.

#### **3. Results**

#### *3.1. Six Clusters of ROIs for the Static FC*

The group adjacent matrix for the static FC is shown in Figure 3a. The horizontal and vertical coordinates represent the sequential numbers of the ROIs. The sequential number and the name of each ROI can be found in Table A1 in Appendix A.

**Figure 3.** (**a**) The group adjacent matrix for the static FC. (**b**) The reorganization of the group adjacent matrix based on the 6 clusters obtained by applying the K-means clustering algorithm to the group adjacent matrix. Since the *i*th row and the *i*th column of the group adjacent matrix are reorganized simultaneously, the reorganized matrix is also symmetric. (**c**) The 6 clusters are shown on a surface rendering of the brain. C1: cluster 1; C2: cluster 2; C3: cluster 3; C4: cluster 4; C5: cluster 5; C6: cluster 6.

Rows of the group adjacent matrix can be clustered into six clusters by the K-means clustering algorithm (Figure 3b). The numbers of rows in clusters 1–6 are 26, 29, 23, 35, 30, and 17, respectively (Table 1). The ROIs in clusters 1–6 can be found in the third and fourth columns of Tables A2–A7 in Appendix A. Since each row of the adjacent matrix corresponds to a ROI, the six clusters can also be shown on a surface rendering of the brain (Figure 3c), which resembles Figure 1a to a large extent.

The average of the squared Euclidean distances from all ROIs in each of the six clusters to the centroid of cluster *i*(*i* = 1, 2, 3, 4, 5, 6) is also evaluated, as shown in Figure 4a–f. For each centroid, among the six averaged distances, the averaged distance from the cluster *i*(*i* = 1, 2, 3, 4, 5, 6) to the centroid of cluster *i* is the lowest. This is consistent with the main idea of the K-means clustering algorithm.

**Figure 4.** The average of the squared Euclidean distances from all ROIs in each of the six clusters to the centroid of cluster *i*(*i* = 1, 2, 3, 4, 5, 6). (**a**) Centroid of cluster 1. (**b**) Centroid of cluster 2. (**c**) Centroid of cluster 3. (**d**) Centroid of cluster 4. (**e**) Centroid of cluster 5. (**f**) Centroid of cluster 6. The error bars represent standard deviations.

#### *3.2. The Overlap Ratios between the Six Clusters for the Static FC and the Six Functional Networks*

The overlap ratios between each cluster and each functional network is analyzed in Table 1. The overlap ratios between cluster 1 and the default network, the frontal-parietal network, the cingulo-opercular network, the sensorimotor network, the occipital network, as well as the cerebellum network are 25/26 (≈96.15%), 0, 1/26 (≈3.85%), 0, 0, and 0, respectively. Obviously, the overlap ratio between cluster 1 and the default network is the highest. Thus, cluster 1 corresponds to the default network. Similarly, we can obtain that clusters 2–6, respectively, correspond to the frontal-parietal network, the cingulo-opercular network, the sensorimotor network, the occipital network, and the cerebellum network, with the overlap ratios, respectively, equaling 20/29 (≈68.97%), 21/23 (≈91.30%), 32/35 (≈91.43%), 22/30 (≈73.33%), and 14/17 (≈82.35%). These overlap ratios are high.


**Table 1.** The number of ROIs in the overlapping part between each functional network and each cluster obtained from the static FC.

#### *3.3. The Consistency Ratios between the Six Clusters for the Static FC and the Functional Networks*

Based on the data shown in Table 1, the consistency between the cluster results and the functional networks can also be evaluated. The consistency ratio between cluster 1 and the default network is 25/(25 + 9 + 1) (≈71.43%), in which 9 is the number of ROIs in the default network but not in cluster 1, and 1 is the number of ROIs in cluster 1 but not in the default network. Similarly, we can obtain that the consistency ratios between cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network are 20/(20 + 1 + 9) (≈66.67%), 21/(21 + 11 + 2) (≈61.76%), 32/(32 + 1 + 3) (≈88.89%), 22/(22 + 0 + 8) (≈73.33%), and 14/(14 + 4 + 3) (≈66.67%), respectively. These consistency ratios are high.

#### *3.4. Six Clusters of ROIs for the SampEn of Dynamic FC*

The group adjacent matrix for the SampEn of dynamic FC is presented in Figure 5a. The horizontal and vertical coordinates stand for the sequential numbers of the ROIs. The sequential number and the name of each ROI can be found in Table A1 in Appendix A.

**Figure 5.** (**a**) The group adjacent matrix for the SampEn of dynamic FC. (**b**) The reorganization of the group adjacent matrix based on the six clusters obtained by applying the K-means clustering algorithm to the group adjacent matrix. Since the *i*th row and the *i*th column of the group adjacent matrix are reorganized simultaneously, the reorganized matrix is also symmetric. (**c**) The six clusters are shown on a surface rendering of the brain. C1: cluster 1; C2: cluster 2; C3: cluster 3; C4: cluster 4; C5: cluster 5; C6: cluster 6.

Rows of the group adjacent matrix can be divided into six clusters by the K-means clustering algorithm (Figure 5b). The numbers of rows in clusters 1–6 are 30, 23, 27, 33, 27, and 20, respectively (Table 2). The ROIs in clusters 1–6 can be found in the fifth and sixth columns of Tables A2–A7 in Appendix A. The six clusters can also be shown on a surface rendering of the brain (Figure 5c), which resembles Figures 1a and 3c to a large extent.

Furthermore, other values of *K*(*K* = 2, ... , 12) are also tried in the K-means clustering algorithm, and the optimal value of *K* is determined by the elbow criterion of the cluster validity index, which is defined as the ratio of within-cluster distances to between-cluster distances [15,20,27]. The dependence of the cluster validity index on *K* is shown in Figure 6. It is seen that two elbows appear at *K* = 4 and 6 due to the changes of slopes of the trend lines. Thus, the optimal values of *K* are 4 and 6. In order to compare the cluster results with the six functional networks already discussed in the literature [25], *K* is fixed at 6 in the present paper.

The average of the squared Euclidean distances from all ROIs in each of the six clusters to the centroid of cluster *i*(*i* = 1, 2, 3, 4, 5, 6) is calculated, as shown in Figure 7a–f. For each centroid, among the six averaged distances, the averaged distance from the cluster *i*(*i* = 1, 2, 3, 4, 5, 6) to the centroid of cluster *i* is the lowest. This is also in line with the main idea of the K-means clustering algorithm.

#### *3.5. The Overlap Ratios between the Six Clusters for the SampEn of Dynamic FC and the Six Functional Networks*

The overlap ratio between each cluster and each functional network is analyzed in Table 2. By evaluating the overlap ratio between each cluster and each functional network, we find that clusters 1–6, respectively, correspond to the default network, the frontal-parietal network, the cingulo-opercular network, the sensorimotor network, the occipital network, and the cerebellum network, with the overlap ratios, respectively, equaling 29/30 (≈96.67%), 20/23 (≈86.96%), 23/27 (≈85.19%), 30/33 (≈90.91%), 22/27 (≈81.48%), and 18/20 (≈90.00%). These overlap ratios are very high.

**Figure 6.** The dependence of the cluster validity index on *K*. The thin solid, dotted, and bold solid lines are trend lines of the filled circles. Since slopes of the trend lines change significantly at *K* = 4 and 6, based on the elbow criterion, the optimal values of *K* are 4 and 6.

**Figure 7.** The average of the squared Euclidean distances from all ROIs in each of the six clusters to the centroid of cluster *i*(*i* = 1, 2, 3, 4, 5, 6). (**a**) Centroid of cluster 1. (**b**) Centroid of cluster 2. (**c**) Centroid of cluster 3. (**d**) Centroid of cluster 4. (**e**) Centroid of cluster 5. (**f**) Centroid of cluster 6. The error bars represent standard deviations.

#### *3.6. The Consistency Ratios between the Six Clusters for the SampEn of Dynamic FC and the Six Functional Networks*

Based on the data shown in Table 2, the consistency ratios between the six clusters obtained from the SampEn of dynamic FC and the six functional networks are evaluated. The consistency ratios between cluster 1 and the default network, cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network are 29/(29 + 5 + 1) (≈82.86%), 20/(20 + 1 + 3) (≈83.33%), 23/(23 + 9 + 4) (≈63.89%), 30/(30 + 3 + 3) (≈83.33%), 22/(22 + 0 + 5) (≈81.48%), and 18/(18 + 0 + 2) (≈90.00%), respectively. These consistency ratios are very high.


**Table 2.** The number of ROIs in the overlapping part between each functional network and each cluster obtained from the SampEn of dynamic FC.

#### *3.7. The SampEn of Dynamic FC is More Effective Than the Static FC*

For the two different measurements (the static FC and the SampEn of dynamic FC), the overlap ratios between cluster 1 and the default network, cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network are shown in Figure 8. For cluster 3, the overlap ratio corresponding to the static FC (91.30%) is larger than that corresponding to the SampEn of dynamic FC (85.19%). For cluster 4, the overlap ratio corresponding to the static FC (91.43%) is slightly larger than that corresponding to the SampEn of dynamic FC (90.91%). For the other four clusters (clusters 1, 2, 5, and 6), the overlap ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC. For clusters 1, 2, 5, and 6, the overlap ratios corresponding to the SampEn of dynamic FC are 96.67%, 86.96%, 81.48%, and 90.00%, whereas the overlap ratios corresponding to the static FC are 96.15%, 68.97%, 73.33%, and 82.35%.

**Figure 8.** The overlap ratios between cluster 1 and the default network, cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network for the two different measurements.

For the two different measurements, the consistency ratios between cluster 1 and the default network, cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network are shown in Figure 9. For cluster 4, the consistency ratio corresponding to the static FC (88.89%) is larger than that corresponding to the SampEn of dynamic FC (83.33%). For the other five clusters, the consistency ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC. For clusters 1, 2, 3, 5, and 6, the consistency ratios corresponding to the SampEn of dynamic FC are 82.86%, 83.33%, 63.89%, 81.48%, and 90.00%, whereas the consistency ratios corresponding to the static FC are 71.43%, 66.67%, 61.76%, 73.33%, and 66.67%.

**Figure 9.** The consistency ratios between cluster 1 and the default network, cluster 2 and the frontal-parietal network, cluster 3 and the cingulo-opercular network, cluster 4 and the sensorimotor network, cluster 5 and the occipital network, and cluster 6 and the cerebellum network for the two different measurements.

According to the results shown in Figures 8 and 9, we conclude that the SampEn of dynamic FC is more effective than the static FC in clustering different ROIs into different functional networks. This phenomenon can be interpreted by evaluating the similarity between the adjacent matrix generated based on the six functional networks (Figure 1b) and the group adjacent matrix for the static FC (Figure 3a) or for the SampEn of dynamic FC (Figure 5a). The similarity is evaluated by the squared Euclidean distance, and a smaller distance means more similarity. The distances from the adjacent matrix shown in Figure 1b to the group adjacent matrices shown in Figure 3a and in Figure 5a are 2409.58 and 2376.52, respectively. The latter is smaller than the former, i.e., the similarity between the adjacent matrix shown in Figure 1b and the group adjacent matrix shown in Figure 5a is larger than the similarity between the adjacent matrix shown in Figure 1b and the group adjacent matrix shown in Figure 3a. This causes the SampEn of dynamic FC to be more effective than the static FC in clustering different ROIs into different functional networks.

#### **4. Conclusions and Discussion**

Different brain regions in the human brain functionally interact with each other to construct multiple functional networks. Identifying the function of each functional network and the brain regions contained in each functional network is very important for understanding the human brain. The present study tests the feasibility of using the K-means clustering algorithm to identify the functional networks based on the FC, including the static FC and the dynamic FC. By applying the K-means clustering algorithm to the static FC or the SampEn of dynamic FC between different ROIs determined by the Dosenbach's template, we show that the Dosenbach's 160 ROIs can be divided into six clusters which show high overlap and consistency ratios with the six functional networks identified by applying modularity optimization on the average FC matrix across a large cohort of healthy subjects. The results indicate that the combination of the K-means clustering algorithm and the FC can identify the functional networks of the human brain. The K-means algorithm has been commonly used to parcellate cortical or subcortical regions based on the static FC [30–38]. These previous studies along with the present study extend the application of machine learning methods in brain sciences.

Furthermore, we show that, for four of six clusters, the overlap ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC, and for five of six clusters, the consistency ratios corresponding to the SampEn of dynamic FC are larger than that corresponding to the static FC. This indicates that nonlinear dynamic characteristics of the FC is more effective than the static characteristics of the FC in identifying brain functional networks. In our previous studies, by characterizing the nonlinear characteristics of dynamic FC in healthy subjects and patients with schizophrenia, we have shown that SampEn of the amygdala-cortical FC in healthy subjects decreased

with age increasing, and the visual cortex of the patients with schizophrenia exhibited significantly higher SampEn than that of the healthy subjects [24,43]. In the future, nonlinear characteristics of dynamic FC should be deeply used to characterize properties of brain functional networks and the complexity of the human brain.

**Author Contributions:** Conceptualization, Y.J. and H.G.; Data curation, Y.J.; Formal analysis, Y.J. and H.G.; Funding acquisition, Y.J. and H.G.; Investigation, Y.J. and H.G.; Methodology, Y.J. and H.G.; Resources, Y.J. and H.G.; Software, Y.J.; Supervision, H.G.; Visualization, Y.J.; Writing-original draft, Y.J.; Writing-review & editing, H.G.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos: 11802086, 11872276, and 11572225) and the Scientific and Technological Project of Henan Province (Grant No: 192102210263).

**Acknowledgments:** We would like to thank www.nitrc.org for allowing us to access the database used in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**


**Table A1.** The names and the sequential numbers of 160 ROIs.


**Table A2.** ROIs in the default network and in cluster 1 for the static FC and the SampEn of dynamic FC. ROIs in cluster 1 but not in the default network are marked by underlines.

**Table A3.** ROIs in the frontal-parietal network and in cluster 2 for the static FC and the SampEn of dynamic FC. ROIs in cluster 2 but not in the frontal-parietal network are marked by underlines.



**Table A4.** ROIs in the cingulo-percular network and in cluster 3 for the static FC and the SampEn of dynamic FC. ROIs in cluster 3 but not in the cingulo-percular network are marked by underlines.

**Table A5.** ROIs in the sensorimotor network and in cluster 4 for the static FC and the SampEn of dynamic FC. ROIs in cluster 4 but not in the sensorimotor network are marked by underlines.



**Table A6.** ROIs in the occipital network and in cluster 5 for the static FC and the SampEn of dynamic FC. ROIs in cluster 5 but not in the occipital network are marked by underlines.

**Table A7.** ROIs in the cerebellum network and in cluster 6 for the static FC and the SampEn of dynamic FC. ROIs in cluster 6 but not in the cerebellum network are marked by underlines.


#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **fNIRS Complexity Analysis for the Assessment of Motor Imagery and Mental Arithmetic Tasks**

#### **Ameer Ghouse 1,2,\*, Mimma Nardelli 1,2 and Gaetano Valenza 1,2**


Received: 10 June 2020; Accepted: 8 July 2020; Published: 11 July 2020

**Abstract:** Conventional methods for analyzing functional near-infrared spectroscopy (fNIRS) signals primarily focus on characterizing linear dynamics of the underlying metabolic processes. Nevertheless, linear analysis may underrepresent the true physiological processes that fully characterizes the complex and nonlinear metabolic activity sustaining brain function. Although there have been recent attempts to characterize nonlinearities in fNIRS signals in various experimental protocols, to our knowledge there has yet to be a study that evaluates the utility of complex characterizations of fNIRS in comparison to standard methods, such as the mean value of hemoglobin. Thus, the aim of this study was to investigate the entropy of hemoglobin concentration time series obtained from fNIRS signals and perform a comparitive analysis with standard mean hemoglobin analysis of functional activation. Publicly available data from 29 subjects performing motor imagery and mental arithmetics tasks were exploited for the purpose of this study. The experimental results show that entropy analysis on fNIRS signals may potentially uncover meaningful activation areas that enrich and complement the set identified through a traditional linear analysis.

**Keywords:** fNIRS; entropy; complexity analysis; nonlinear analysis; brain dynamics; mental arithmetics; motor imagery

#### **1. Introduction**

Functional near-infrared spectroscopy (fNIRS) is a noninvasive technique that has found success in analyzing brain function through the lens of metabolic processes and neurovascular coupling [1,2]. Common methods found in the literature analyze fNIRS signals with the assumption that an underlying linear system generated their time series [3]. Though these approaches may find success in some domains, linearity is an ideal assumption when investigating brain physiology. In fact, many physiological systems exhibit nonlinear behavior, meaning there can be further interaction between variables in a system beyond a superposition effect while also having dynamics that the system sub-components may not show. Beyond nonlinearity, physiological systems may exhibit complex dynamics as a result of feedback loops that arise from homeostasis regulation with consequent extreme sensitivity to the system state condition [4–6].

Prior literature has shown that nonlinearities are particularly present in the brain and its related metabolic processes. Functional magnetic resonance imaging (fMRI) and fNIRS data were demonstrated to follow a nonlinear saturating impulse response model [7], and physiological models of cerebral blood

flow dynamics include complex feedback loops between ion channels, metabolism, energy demand, and oxygenation [8]. Furthermore, dynamics of the intrinsic parameters, such as the electrophysiological process that drives neurovascular coupling, also exhibit nonlinear and complex behavior [9,10].

Such nonlinearities found in metabolic processes imply that standard linear models and metrics quantifying linear dynamics defined in the time and frequency domains may potentially underrepresent the physiological processes sustaining functional activity. To this end, entropy can be a powerful tool to characterize a system's regularity or complexity [11]. When applied to the topology of attractors describing a dynamical system in phase space, entropy leads to a robust estimation of regularity of state space evolution, also known as the Kolmogorov–Sinai metric [12]. By exploiting Takens' theorem and the concept of characterizing an attractor through its topological entropy, several algorithms have been developed to find a value that converges to the Kolmogorov–Sinai entropy metric for regularity. Such algorithms include sample entropy (SampEn) [13] and fuzzy entropy (FuzzyEn) [14], which are able to characterize a system's regularity at a single time scale level [15]. On the other hand, metrics, such as distribution entropy (DistEn) [16], have been shown to provide complexity estimates of the system under study.

While entropy analysis has been a widely investigated tool for studying electrophysiological signals, there is a dearth of studies regarding entropy applied to metabolic processes, as observed in fNIRS signals. Permutation entropy, i.e., entropy of a time series from an ordinal transform on the continuous data [17,18], has been exploited by Gu et al. to investigate the complexity of fNIRS signals in children affected by attention deficit disorder during working memory tasks [19]. Furthermore, Jin et al. investigated permutation entropy to analyze differences in experts and novices solving science problems [20]. Studying frontal cortex fNIRS signals, SampEn was suggested as a biomarker for Alzheimer's disease diagnosis [21–23], and Angsuwatanakul et al. [24] investigated the effects of working memory experiments on SampEn estimated from fNIRS series. Also, though applied as an information theoretic approach to investigate linear effects in fNIRS rather than analyze topological entropy in phase space, differential entropy has been investigated in Keshmeri et al. as a biomarker that preserves variational information in the assessment of working memory [25,26].

Although there is literature for entropy applied to fNIRS signals, there has yet to be an analysis of its regularity and complexity during standard cognitive load tests, such as motor imagery and mental arithmetics. Besides, previous studies using entropy were not performed using a time stamped controlled block design protocol. Thus, it is not yet clear how well entropy as an estimate works when activity is controlled in time. Furthermore, a comparison with standard methods deserves scrutiny. To overcome these oversights, this study aims to uncover SampEn, FuzzyEn, and DistEn estimates of hemoglobin, deoxyhemoglobin, and total hemoglobin in mental arithmetics and motor imagery experiments in order to perform a comparison with traditional methods in fNIRS signals analysis. Concretely, we hypothesize that by considering nonlinear and complex characterizations of metabolic processed observed in fNIRS signals, more information, as expressed by cortical activity correlates, can be gleaned regarding physiological and psychophysiological phenomena than what can be considered using only linear analyses. For the purpose of this study, we used an open access dataset provided by Shin et al. [27], whose details on methodology and results follow below.

#### **2. Materials and Methods**

#### *2.1. Block Design*

The dataset used in this study is openly available and fully described in [27]. Briefly, twenty-nine subjects (aged 28.5 ± 3.7, 15 females) were involved in the experiment. Left and right hand motor imagery constituted one set of trials performed, and the other set of trials were baseline and mental

arithmetics. There were three trials of each of the aforementioned experiments per subject. fNIRS and electroencephalography (EEG) series were acquired simultaneously during the whole duration of the experiment using 30 EEG channels and 36 fNIRS channels, and the sampling rate for fNIRS signals was 10 Hz. The 36 fNIRS channels were resolved from a set of 14 sources to 16 detectors matching as illustrated in Figure 1.

The experimental protocol began with a 60 s rest, after which subjects were presented an instruction (either a "←" or "→" for motor imagery experiments, and either a "-" or an arithmetic task in comparing baseline vs mental arithmetic) on the screen telling them which task were to be performed. Afterwards, the individual performed the task for 10 s, with a subsequent 15 s rest before the next task. After 20 repetitions of these instructions and tasks, a 60 s ending rest was performed. Mental arithmetic/baseline trials were performed independently from motor imagery trials.

**Figure 1.** Position of the Optodes. Positions labeled with "D" refer to detectors while positions labeled with "S" are sources. The lines demonstrate coupling between sources and detectors.

#### *2.2. Hemoglobin Extraction from fNIRS Signals*

In continuous wave fNIRS acquisitions, light radiations from two different wavelengths are used to create a system of equations that can resolve hemoglobin content. These wavelengths are generally chosen to be in the range of the physiological window where water and hemoglobin absorption is particularly low (650 nm to 1350 nm). To this extent, the "modified Beer–Lambert law" provides a mathematical

expression relating absorption measured with a detector and the concentration of a chromophore as seen in Equation (1) [28]:

$$\mu\_a(\lambda, t) = \log(\frac{I\_o(\lambda)}{I(\lambda, t)}) = \sum\_{i=0}^n c\_i(t)\varepsilon\_i(\lambda)\rho DPF + G \tag{1}$$

where *μ<sup>a</sup>* is the absorption coefficient at a given wavelength *λ* and time *t*, *Io* is incident light intensity, *I* is the detected intensity that changes with time, *c* are the chromophore concentrations of interest, *ρ* is the separation between a light source and detector, DPF is the correction factor for a best estimate of a light path through a tissue, and G is the loss of light due to scattering. In a continuous wave setting, differential concentrations Δ*c*, related by differential absorption Δ*μa*, are the parameters that are analyzed in the fNIRS signals. This allows for a significant simplification of the expression above when assuming that scattering loss is a constant in time, yielding differential concentrations that can be resolved through a simple linear system of equations given multiple wavelengths, as seen in Equation (2):

$$\Delta\mu\_d(\lambda, t) = \log(\frac{I\_b(\lambda)}{I(\lambda, t)}) = \sum\_{i=0}^{n} \Delta c\_i(t)\varepsilon\_i(\lambda)\rho DPF \tag{2}$$

where *Ib* is the intensity detected at a baseline of interest.

From the modified Beer–Lambert law, the differential concentration of deoxyhemoglobin can be retrieved by choosing two wavelengths on opposing sides of the isobestic point of the absorption spectra of oxyhemoglobin and deoxyhemoglobin and solving a linear system of equations. In the methods of Shin et al., 760 nm and 850 nm were used as wavelengths. In this study, an evaluation on such differential concentrations as well as total Hb during activity will be performed.

#### *2.3. fNIRS Data Preprocessing*

Figure 2 shows an overview of the preprocessing pipeline. The signal was transformed from optical densities into Hb and HbO using the modified Beer–Lambert law. For the modified Beer–Lambert law, the first 60 s were considered as a baseline, corresponding to the resting state. A first Butterworth lowpass filter with a cutoff frequency at 0.6 Hz and filter order 6 was applied to fNIRS data to highlight the hemodynamic response. This was considered as 0.6 Hz is the upper end of cut-off frequencies used in literature [29]. This is significant for preserving the full dynamics of hemoglobin, including the high frequency components, which can uniquely affect the topology of the attractor in phase space and render different estimates of entropy. In hand, we must accept the risk of physiological phenomena, such as Mayer waves, contaminating the entropy estimates. A second band-pass filter with cutoffs 0.8 Hz and 2 Hz and filter order 6 was used to capture pulsatile dynamics of hemoglobin [30]. Afterwards, a wavelet filtering approach was used to further reduce noise, particularly related to motion artifacts, in the oxy- and deoxyhemoglobin signals [31]. This wavelet filtering approach works by decomposing the time series into nine levels using a daubechies five mother wavelet, subsequently thresholding detail coefficients that have low probability (*p* < 0.1) given the detail coefficients are sampled from a normal distribution. After the wavelet filter, the time series were separated into epochs representing blocks of activity. Each channel at each activity block was differentially referenced to the mean of the previous 5 s of said channel. The data was then further processed to extract features such as entropy and mean values of hemoglobin.

**Figure 2.** Pipeline for processing functional near-infrared spectroscopy (fNIRS) data.

#### *2.4. Entropy Analysis*

The entropy metrics SampEn, FuzzyEn, and DistEn were extracted as regularity and complexity characterizations of fNIRS data. For each fNIRS signal (Hb, HbO, and total hemoglobin) and for the multivariate embedding derived from a concatenation of the Hb and HbO embedding, the optimal time delay was chosen as the first zero of the autocorrelation while the optimal embedding dimension was found using the false nearest neighbours algorithm [32].

To create the embeddings, we started from a time series x(t) of *N* samples. Having determined the time lag *τ* and embedding dimension *m*, the states *X<sup>m</sup> <sup>t</sup>* of the reconstructed attractor can be represented in vector form as follows:

$$X\_t^{\text{wt}} = \{ \mathbf{x}(t), \mathbf{x}(t+\tau), \dots, \mathbf{x}(t+(m-1)\tau) \}\tag{3}$$

When reconstructing an attractor using several variables, i.e., the concatenated attractor (Concat), the above expression is modified in the following way:

$$X\_t^m = \{ \mathbf{x}(t), \mathbf{x}(t+\tau), \dots \mathbf{x}(t+(m-1)\tau), y(t), y(t+\tau), \dots, y(t+(m-1)\tau) \}\tag{4}$$

From the reconstructed attractor, entropy estimates may be computed. For SampEn, the estimate can be obtained, as follows:

$$\text{SampEn} = -\log(\frac{\sum\_{i=1, i \neq j}^{N-m} \frac{1}{N-m-1} \text{number of } |X\_i^{m+1} - X\_j^{m+1}| < R}{\sum\_{i=1, i \neq j}^{N-m} \frac{1}{N-m-1} \text{number of } |X\_i^m - X\_j^m| < R}) \tag{5}$$

where *R* refers to a user set deviance of states to binarize the distance metric. In this study, it is set to 0.2 ∗ *σx*, a setting widely used in previous studies with theoretical justifications, where *σ<sup>x</sup>* is the standard deviation of the considered time series [33,34].

FuzzyEn uses a fuzzy membership function instead of the heaviside function to calculate the correlation integral. In this study, we employed an exponential function, as follows:

$$\boldsymbol{\phi}^{\rm m} = \frac{1}{N - m} \sum\_{i=1}^{N-m} \frac{1}{N - m - 1} \sum\_{j=1, j \neq i}^{N-m} \boldsymbol{\varepsilon} \frac{(-|\boldsymbol{\lambda}\_i^m - \boldsymbol{\lambda}\_j^m|)^K}{\boldsymbol{\kappa}} \tag{6}$$

The value of *K* was set to 2. Afterwards, FuzzyEn can be derived as the ratio between the above fuzzy function with the result of a fuzzy function of an order greater [35].

$$FuxzyEn = -\log(\frac{\phi^{m+1}}{\phi^m})\tag{7}$$

DistEn is less dependent on parameter selection in comparison to FuzzyEn and SampEn, given that the parameter *R* is no longer required. A histogram is constructed from the distance matrix, and the Shannon entropy of the empirical probability density function is computed. To make the algorithm faster, we extracted the upper triangle of the distance matrix, as it should be symmetrical, meaning that lower triangle contains redundant information. Additionally, the diagonal is removed from the entropy estimate as it should be a zero vector when considering that the self-similar distance is zero. Bin size was estimated by using Scott's method [36].

#### *2.5. Statistical Analysis*

A bootstrapped third moment test was performed with linear time series surrogate samples generated by an amplitude adjusted Fourier transform of the original time series and phase scrambling in order to test the null hypothesis that the original time series was generated from a linear system [37]. Two-hundred surrogate series were generated in order to determine a p-value, with the third moment calculated for *τ* = 1 lag as illustrated in the following equation [38].

$$t^{c3}(\tau) = \preccurlyeq \mathbf{x}\_k \cdot \mathbf{x}\_{k+\tau} \cdot \mathbf{x}\_{k+2\tau} \text{ } \tag{8}$$

The percentage of significant time series for each channel and each parameter, either Hb, HbO, or total Hb, are shown in the results. Significance is determined using an *α* = 0.05.

After ascertaining nonlinearity, further non-parametric tests were performed on entropy and mean hemoglobin results when considering the non-Gaussian distribution of the metrics. Friedman non-parametric statistical tests for paired data were performed in order to determine whether repetitions of activities in each trial were significantly different. Afterwards, a Friedman test was applied using a median summary statistic over trials to compare significant cortical areas of activation between the four tasks (i.e., baseline, mental arithmetic, left hand, and right hand motor imagery). Multiple comparison tests were then performed between pairs of tasks using Wilcoxon signed rank tests for paired data, and the statistical significance was set to 0.05 when considering a Bonferroni correction rule over the four different activity comparisons.

Group-wise and channel-wise multiple comparison results for each metric are displayed using both p-value topographic maps and topographic maps displaying Δ value differences between tasks for a given metric. Cortical regions in the topographic maps that are not covered by the optodes, as seen in Figure 1, are inferred using a bilinear interpolation.

#### **3. Results**

#### *3.1. Nonlinearity Test*

As illustrated in Figure 3, an analysis of the third moment for each time series in Hb, HbO, and total hemoglobin demonstrates that the majority of time series exhibits nonlinear behavior, rejecting the null hypothesis that a linear system generated the time series.

**Figure 3.** Topographic maps from channel-wise third moment tests displaying the fraction of time series from each channel having statistical significance, where the colorbar indicates the value of the fraction.

#### *3.2. Analysis of Repetitions within Tasks*

Through the Friedman statistical test on repetitions, it can be seen by Table 1 that we were able to accept the null hypothesis that there were no significant differences between the repetitions for either mental arithmetic, left hand imagery, right hand imagery, or baseline when using any of the statistics of mean, SampEn, or DistEn over any set of hemoglobin time series representation. On the other hand, we could reject the null hypothesis for the FuzzyEn comparisons in the case of using total hemoglobin time series and the multivariate topology reconstructed from both oxyhemoglobin and deoxyhemoglobin.

**Table 1.** Table of statistical power *p*-values from the Friedman analysis. *p*-values are bonferroni corrected. \* denotes that using an alpha of 0.01 we must reject the null hypothesis that there were no significant variations between repetitions. This particularly occurs for FuzzyEn in the total and the concatenated case for deoxyhemoglobin.


Given this result, subsequent post-hoc analyses focused on mean estimates, SampEn, and DistEn for each time series. Furthermore, when considering that the repetitions of these metrics did not show significant differences, the median value of each estimate over repetitions was used as a summary statistics for further inter-subject analyses.

#### *3.3. Between-Task Statistical Analysis*

Cortical areas with significant statistical differences between baseline, mental arithmetic, right hand, and left hand motor imagery tasks according to a Friedman test analysis on mean, SampEn and DistEn analyses can be seen in Figure 4. Estimates on the oxyhemoglobin signal showed overlapping areas of significance between mean estimate and both entropy estimates in the occipital regions. On the other hand, DistEn estimates on deoxyhemoglobin signal had significant changes between tasks over the somatosensory cortex that were not exhibited in the mean estimate. For the total hemoglobin signal, both SampEn and DistEn unraveled further information that mean estimates could not, where DistEn exhibited significant changes in the occipital area, and SampEn exhibited changes in the parietal area. In the concatenated topology, DistEn and SampEn exhibited different subsets of cortical activations.

**Figure 4.** *p*-value topographic maps from channel-wise Friedman tests displaying significant statistical differences between all tasks in the experimental protocol (baseline, mental arithmetic, right hand, and left hand motor imagery). Y (green) areas indicate where we could reject the null hypothesis that activity was the same in all the tasks, whereas N (white) areas indicate where we could not reject the null hypothesis.

#### *3.4. Multiple Comparison Analysis*

Figure 5 shows cortical areas that were associated with significant statistical differences between mental arithmetic activity and baseline activity for a given estimate according to Wilcoxon non-parametric tests. When analyzing oxyhemoglobin, SampEn displayed regions in the occipital cortex that were not highlighted by the mean estimates; DistEn did not seem to add further information. From deoxyhemoglobin signal analysis, DistEn uncovered significant changes over the left occipital region that were unobserved in the mean estimate analysis; SampEn did not seem to add new information. On total hemoglobin, significant changes between tasks were found in the right occipital cortex from DistEn, which were unobserved in mean estimates. From the concatenated signal, DistEn displayed information in the parietal cortex that was unobserved by previous analysis. From visual inspection on Figure 5, it seemed that mental arithmetic activity was generally associated with higher mean and a lower irregularity and complexity levels than baseline.

*Entropy* **2020**, *22*, 761

**Figure 5.** *p*-value topographic maps from channel-wise Wilcoxon non-parametric tests displaying significant statistical differences between mental arithmetic activity and baseline activity. Y (green) areas indicate statistically significant changes between tasks, whereas N indicates non-significant changes. The colormap topoplots display estimate differences between baseline (B) and mental arithmetic (M) tasks, with red indicating higher values for mental arithmetic than baseline and blue indicating lower values for mental arithmetic as compared to baseline.

Further Wilcoxon tests were performed to show cortical areas that were associated with a significant statistical difference between left hand motor imagery and baseline for a given estimate, as seen in Figure 6. When analyzing oxyhemoglobin, both DistEn and SampEn showed significant changes between tasks over a larger region than the mean estimate, especially in the right occipital cortex. Furthermore, deoxyhemoglobin activity in the left temporal and sensorimotor cortices was highlighted by both entropies. A total hemoglobin analysis confirmed that DistEn and SampEn highlight further changes that were not seen in a mean estimate analysis. With visual inspection on Figure 6, it appeared that SampEn inversely mapped mean estimate changes over the the frontal, motor, and parietal regions for the oxyhemoglobin signals. For deoxyhemoglobin, higher mean estimates over the right hemisphere were associated with left hand motor imagery activity. SampEn increased over the frontal areas during left hand motor imagery tasks with respect to baseline with no changes over the posterior areas. Changes in total hemoglobin signal seemed similar to oxyhemoglobin.

*Entropy* **2020**, *22*, 761

**Figure 6.** *p*-value topographic maps from channel-wise Wilcoxon non-parametric tests displaying significant statistical differences between left hand imagery activity and baseline activity. Y (green) areas indicate statistically significant changes between tasks, whereas N indicates non-significant changes. The colormap topoplots display estimate differences between baseline (B) and left hand imagery (L) tasks, with red indicating higher values for left hand imagery vs baseline and blue indicating lower values for left hand imagery vs baseline.

From Figure 7, another set of Wilcoxon non-parameteric test results can be seen, showing cortical areas that were associated with a significant statistical difference between right hand motor imagery and baseline for a given estimate. While mean estimates were associated with few significant changes between tasks, SampEn and DistEn showed significant differences over several areas, especially in a oxyhemoglobin and total hemoglobin analysis. Particularly, in a oxyhemoglobin analysis, DistEn showed significant changes over the frontal, right, and left occipital areas, which were complemented by further changes over parietal cortices by SampEn. For deoxyhemoglobin signal, complementary parietal activity appeared in DistEn while SampEn changes were a subset of the mean estimates. In the case of total hemoglobin, changes over the sensorimotor and parietal cortices were found using SampEn, while DistEn and mean estimates did not show significant changes between tasks.

Using further visual inspection analysis on Figure 7, the trend appears to be that higher mean, SampEn, and DistEn values over the frontal areas were more associated with right hand motor imagery activity, whereas higher estimate values over the posterior areas were associated with baseline activity. In the case of deoxyhemoglobin, higher mean estimates over the right hemisphere were associated with

right hand motor imagery activity. SampEn increased over the frontal areas during right hand motor imagery tasks with respect to baseline with no changes over the central posterior areas. Changes in total hemoglobin signal seemed similar to oxyhemoglobin ones.

**Figure 7.** *p*-value topographic maps from channel-wise Wilcoxon non-parametric tests displaying significant statistical differences between right hand imagery activity and baseline activity. Y (green) areas indicate statistically significant changes between tasks, whereas N indicates non-significant changes. The colormap topoplots display estimate differences between baseline (B) and right hand imagery (R) tasks, with red indicating higher values for right hand imagery vs baseline and blue indicating lower values for right hand imagery vs baseline.

Figure 8 shows cortical areas that were associated with a significant statistical difference between right hand and left hand motor imagery for a given estimate according to Wilcoxon non-parametric tests. Complementary left occipital activity was uncovered by DistEn for oxyhemoglobin, while a deoxyhemoglobin analysis using SampEn uncovered unique parietal activity changes between tasks. In the case of total hemoglobin, larger parietal changes were found in DistEn than the mean estimate, while SampEn exhibited changes in the right temporal regions. Visual inspection analysis on Figure 8 shows a trend of left hand motor imagery activity being associated with higher mean, irregularity and complexity levels than right hand motor imagery activity over the frontal areas, while an opposite trend seemed to be observed over the posterior regions. Particularly, changes over the frontal cortex in mean estimates seemed similar to SampEn differences in oxyhemoglobin, while they appeared to be

inversely distributed in DistEn. In deoxyhemoglobin, no differences between left and right hand motor images seemed to occur over the posterior regions in SampEn, whereas DistEn appeared to show similar differences as mean estimates.

**Figure 8.** *p*-value topographic maps from channel-wise Wilcoxon non-parametric tests displaying significant statistical differences between left hand imagery and right hand imagery activities. Y (green) areas indicate statistically significant changes between tasks, whereas N indicates non-significant changes. The colormap topoplots display estimate differences between left hand imagery (L) and right hand imagery (R) tasks, with red indicating higher values for right hand imagery than left hand imagery and blule indicating lower values for right hand imagery than left hand imagery.

#### **4. Discussion**

We investigated changes in fNIRS entropy during mental arithmetics and motor imagery tasks and compared the results with fNIRS standard analysis metrics. Our aim was to test whether entropy analysis could unravel changes in cortical areas that may not be highlighted while using traditional methods that analyze the signal in the time domain. Particularly, we assessed statistical differences in fNIRS signal entropy in four different tasks (baseline, mental arithmetic, right hand, and left hand motor imagery), and compared different entropy metrics—specifically SampEn, FuzzyEn, and DistEn—together with mean value estimates of hemoglobin, deoxyhemoglobin, and total hemoglobin.

Previous studies used entropy estimates in protocols of long time windows with unspecified timing of events in the signal, as in the case of cognitive capacity analysis in Alzheimer's [21–23]. Nevertheless,

regularity and complexity analyses of fNIRS signals during standard cognitive load tests, such as motor imagery and mental arithmetics, were not investigated to the best of our knowledge.

Through a test of nonlinearity, we were able to ascertain that the majority of the considered time series demonstrated nonlinear behavior. Nonlinearity testing was necessary for validating whether quantifying the extent of nonlinear behavior could be a value of interest in the analysis of functional activity. To that extent, our analysis corroborated studies performed in the past, such as the evidence demonstrated in Khoa et al., where they performed similar nonlinearity tests [38].

Our study showed that FuzzyEn applied to total hemoglobin and the concatenated attractor from the open dataset demonstrated significant differences between task repetitions (see Table 1); therefore, only SampEn and DistEn were retained for further analyses on fNIRS regularity and complexity at a task level. In fact, this result allowed for subsequent comparison analyses between the four tasks to be performed using a median summary statistic for entropy and the mean estimates over the repetitions rather than considering each repetition independently.

Over the general set of results, complementary areas of functional activity were found in both SampEn and DistEn when compared to mean estimates, as demonstrated in Figures 5–8. For example, in the comparison between mental arithmetics and baseline activities in Figure 5, SampEn was able to uncover particular parietal activity in oxyhemoglobin that mean estimates, using any of the three hemoglobin concentrations (Hb, HbO, THb), were unable to resolve. Furthermore, it appears that both entropy estimates are more sensitive to temporal cortex activity, as seen in Figures 6 and 7, when analyzing motor imagery tasks compared to baseline.

Previous studies highlighted hemodynamic changes during mental arithmetic tasks primarily over the bilateral intraparietal, inferior temporal, and dorsal prefontal sites [39,40]. SampEn and DistEn were both successful in recovering those activity areas as demonstrated in Figure 5. Particularly, SampEn applied to oxyhemoglobin showed changes over parietal structures while deoxyhemoglobin revealed changes over frontal cortical sites. With DistEn, the concatenated series displayed changes over both parietal, frontal, and temporal activity. However, the mean estimates were not able to uncover the parietal cortex changes, but instead were only sensitive to frontal cortex and temporal cortex activity. In the case of mental arithmetics, these results suggest that entropy estimates may be more sensitive to cortical hemodynamic changes than mean estimates given the sample size available. This may be due to the additional quantification of nonlinear and complex dynamics provided by entropy analysis. Where linear effects subside or may not be as significant, nonlinear, and complex behavior may still persist. This could be explained by models that demonstrate short term stimuli resulting in nonlinear behavior in the hemodynamic response [7]. Speculatively, stimuli may become less frequent or last for shorter durations when a subject experiences fatigue from a protocol or has become habituated.

In light of motor imagery tasks, Figure 8 demonstrates that we were able to find activity areas in the expected sensorimotor cortex while using either entropy or mean estimate analysis. Explicitly, both DistEn applied to deoxyhemoglobin and SampEn applied to the concatenated attractor unraveled these expected changes. Furthermore, we observed a lateralization effect in DistEn applied to oxyhemoglobin and the concatenated attractor, as well as SampEn applied to oxyhemoglobin. These results are in accordance with previous findings [41]. This suggests the presence of complementary information supplied by regularity and complexity analysis on fNIRS series. In light of these significant results, it is important to mention that a bilinear spatial interpolation was performed on the topographic maps as mentioned in Section 2.5, thus there could be errors in drawing the true cortical location of activity. It would be important to use a higher density fNIRS cap in future experiments in order to better pinpoint the true cortical location of a specific activity.

The success of entropy estimates in unraveling complementary areas was particularly surprising when considering that the experiments studied here were tailored to leverage strong activations that arise

from a saturating superposition effect, i.e., linear superposition. As mentioned above, it is possible that short-term stimuli were introduced when either the subject became habituated or fatigued. Furthermore, there may also be significant oscillatory behaviors that contribute to the observed nonlinearities observed in the hemodynamic signal that mean value analysis can not detect. For example, changes in pulsatility in the microvessels that arise from cardiac pulses and physical properties of the microvessels may nonlinearly affect the oxygen extraction from the capillaries to the tissue [10].

As has been mentioned in the introduction, biological systems exhibit a vast array of feedback and compensatory loops in order to regulate homeostatic behavior at a neurolobiological level [4–6]. This knowledge brings light to the significance of the study we have presented in leveraging the information in phase space that this complex system projects in the fNIRS time series. However, a clear limitation of the study is that it is purely exploratory, rather than explanatory for the neurobiological activity that underlies the complex system the entropy estimates assess. Nonetheless, this study holds a beacon for future research to investigate the intrinsic complex neurobiological correlates that comprise activity in mental arithmetic and motor imagery tasks.

A natural extension of this study in the future can be to apply fNIRS regularity and complexity analysis to block-free paradigms, such as a clock drawing test [42], or tests that stimulate more complex dynamics related to emotional response [43–45]. Because SampEn was not applied using a multiscale algorithm, future studies can also investigate fNIRS dynamic activity while using a multiscale entropy analysis. Such sophisticated methodology may further highlight complex changes that may be induced by activity on different time scales, such as cardiac pulsatility, arterial blood pressure induced mayer waves, or other nonlinearities driving the hemodynamic response [10,46,47]. Furthermore, in future studies, a dataset using an fNIRS system that includes short source-detector separation channels can be analyzed to regress out artifacts due to skin-blood flow induced changes in the fNIRS signals.

#### **5. Conclusions**

A novel investigation into the analysis of entropy in metabolic processes measured by fNIRS on controlled block design experimental protocols was presented in this study. We conclude that entropy may uncover areas that yield neuronal correlates and that agree with traditional methods of analyzing neuronal correlates while also providing novel complementary areas not seen in mean estimates. Furthermore, entropy estimates seemed to exhibit greater sensitivity with sample size to activity than mean estimates in mental arithmetics. These results shed light on not only the validity, but also the efficacy of using entropy to investigate functional neural activations.

**Author Contributions:** A.G. processed and analyzed the data under the supervision of G.V. and M.N.; Writing of the original draft was performed by A.G.; All authors reviewed and revised the final manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research leading to these results has received partial funding from the European Commission—Horizon 2020 Program under grant agreement n◦ 813234 of the project "RHUMBO" and by the Italian Ministry of Education and Research (MIUR) in the framework of the CrossLab project (Departments of Excellence).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**




c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
