A total of 32 healthy younger people (15 males, 17 females; average age, 23.9 years; standard deviation [SD], 4.7 years; age range, 20–35 years) and 18 healthy older people (7 males, 11 females; average age, 57.5 years; SD, 4.7 years; age range, 51–67 years) were enrolled in this study. These groups were sex-matched (χ² = 0.30, p = 0.59). The older participants were all non-smokers and were not on any medications. Participants with medical or neurological conditions (including epilepsy or head trauma occurring in the past), as well as those with a history of alcohol or drug dependence, were excluded from the current study. All the participants provided their written informed consent following an explanation of study procedures as well as risks and benefits by study personnel. This study was approved by the Ethics Committee of Kanazawa University and was conducted in accordance with the principles of the Declaration of Helsinki and its later amendments. The EEG data used in this study evaluated the dynamics of phase synchronization (Nobukawa et al., 2019a).

Methods for recording and pre-processing EEG data have been reported and established in previous research (Mizuno et al., 2010). Specifically, the participants in the current study sat in a soundproof recording room, and their EEG was measured under controlled room lighting conditions. For EEG measurement, 16 electrodes (Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, Fz, Pz, T5, and T6) were used; this system was based on the recommended electrode arrangement under the international 10–20 system. EEG activity was measured with reference to the binaural connection. The EEG-4518 monitor used for electroencephalogram measurements in this study was manufactured by Nihon Kohden Co., Ltd. (Tokyo, Japan). The sampling frequency was 200 Hz for the recording. The electrode/skin conductance impedance was controlled to within less than 5kΩ for each electrode. Participants' electroencephalogram signals were measured for 10–15 min in a resting state with the eyes closed. Researchers visually inspected the participants' arousal using a video surveillance system; participants were asked to close their eyes, and researchers confirmed that only awake epochs were measured. If the alpha and theta oscillations became weaker or stronger, compared with ones at beginning stage of the recording, this duration was not used for evaluation, because this duration belonged to the light sleep stage. Additionally, the EEG signals were visually assessed to identify artifacts, such as muscle activity, blinks, and eye movement; consequently, 60-s (12,000 data points) of artifact-free time-series segments within the EEG signals recorded in the awake state with eyes closed were identified. For each epoch, bandpass filtering with the range of 2.0–60 Hz was applied. The first and last 5-s period (1,000 data point) in each bandpass-filtered epoch were removed to avoid transient behaviors produced by the bandpass filtering process. MSE and MF analyses were performed for 50 consecutive seconds (i.e., 10,000 data points) of epochs. In the PLI analyses, values decrease with increasing epoch length (Fraschini et al., 2016); therefore, it is difficult to identify changes with an increasing epoch length. In addition, using short epoch lengths makes it impossible to capture information on slow frequency components. To balance these considerations, the PLI analysis divided 50 consecutive seconds (10,000 data points) into 10 epochs of 5 s each (Takahashi et al., 2017, 2018; Nobukawa et al., 2020a,b).

The overview of flow for multi-fractal analysis is shown in Figure 1A. In MF analyses, wavelet leaders derived from the coefficients of the discrete wavelet transform are widely used (Jaffard et al., 2006; Wendt and Abry, 2007). MF analysis is an analysis method that uses the Hölder index to represent the fractal dimension of the partial structure that characterizes the structure of data X via spectrum data. The discrete wavelet coefficient of the discrete signal X(t) is given by Equation (1).

Here, ψ₀ is a compact-supported mother wavelet function. The Equation (2) shows one-dimensional wavelet leaders which are time- or frequency-localized suprema of the absolute value of the discrete wavelet coefficients d_X(j, k):

Here, λ=λj,k=[k2-j,(k+1)2-j] represents the time interval of the scale 2^−j. Additionally, 3λ_j,k−1 = ∪λ_j,k ∪ λ_{j, k+1} represents the adjacent time (Wendt and Abry, 2007). The spectrum of singularity of L_X is defined by Equation (3) with wavelet leaders (Jaffard et al., 2006; Wendt and Abry, 2007).

Here, h indicates the Hölder index. Also, q indicates the moment of the wavelet leaders. The scaling index ζ_L(q) is defined by Equation (4). The wavelet leader structure function S_L(q, j) is defined by Equation (5).

Here, n_j indicates the number of samples of X when the scale is 2^j. As the Hölder index h approaches 1.0, the time-series shape becomes more differentiable. However, as the Hölder index h approaches 0, the time-series shape becomes nearly discontinuous. A signal is monofractal if the scaling index ζ_L(q) is a linear function and D(h) converges to a particular h. Contrastingly, the fact that the signal is multi-fractal indicates a scaling index, where ζ_L(q) deviates from linearity and D(h) is widely distributed in h. In this study, to capture the profile of D(h), the primary cumulant c₁ of D(h) was used as an indicator of the smoothness of the entire time series signal, and the secondary cumulant c₂ was used as an index evaluating the local fluctuation of the time-series signal. For the multi-fractal time series, D(h) is distributed around c₁. Therefore, the degree of distribution of D(h) reflects the multi-fractal property, which corresponds to |c₂|. The time-series with large (small) multi-fractality (|c₂|) exhibits intermittent and transient behavior with large (small) amplitude (Ihlen, 2012); while, the complexity notified by c₁ reflects the degree of complexity for temporal behavior in entire time-range, instead of intermittent behavior (see Supplementary Material). In this study, multi-fractal analysis was performed using the wavelet toolbox in MATLAB (https://jp.mathworks.com/products/wavelet.html; MathWorks, Natick, MA, USA).

The overview of flow for MSE analysis was shown in Figure 1B. MSE analysis was used to assess the temporal scale dependence for EEG time series complexity (Costa et al., 2002). The time-series sample entropy of the random Z-score variables {x₁, x₂, ..., x_N} is given by Equation (6).

Here, C_m(r) is the probability of |xim-xjm|<r(i≠j,i,j=1,2,...) among all pairs of i and j. xim indicates an m-dimensional vector xim={xi,xi+1,...,xi+m-1}. In the MSE analysis, {x₁, x₂, ..., x_N} is calculated using Equation (7) for coarse-grained time series y_j.

Here, {y₁, y₂, ..., y_N} represents the observed signals. τ(τ = 1, 2, ...) represents the temporal scale. In this study, we set m = 2 and r = 0.2 (Costa et al., 2002) and MSE analysis was performed using the Physio Toolkit toolbox in MATLAB (http://physionet.incor.usp.br/physiotools/sampen/).

The PLI was obtained to measure phase synchronization, and the characteristics of the synchronization signal were quantitatively estimated. The EEG signal was divided into five frequency bands: the delta (2–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma bands (30–60 Hz). Here, several studies showed that artifacts in the gamma band due to muscle activity is larger compared with slower frequency bands (Whitham et al., 2007, 2008). Therefore, through visual examination of EEG signals, the time including muscle activity was avoided in the evaluation epochs. Each band division divides the signal at time t, and the point ϕ_a is indicated by phase ϕ_a(t) and amplitude Aa(t) using the Hilbert transform. In addition, the phase difference Δϕ_ab(t_i) observed between signals with two different points a and b at time t_i is given by Equations (8) and (9) (Stam et al., 2007).

The PLI of the signal between the two points a and b for the duration T is given by Equation (10).

When signals with the same source are observed at different points, Δϕ_ab(t_i) is 0 and Δϕ_mod(t_i) = 0; subsequently the PLI_ab value becomes 0. In addition, the observation at the point on the opposite side of the electric dipole is defined as Δϕ_ab(t_i) = π within Equation (8) in cases where the signal source is assumed to follow the dipole model. This yields PLI_ab = 0. The average PLI of any electrode a via another electrode b = 1, 2, ..., K(b ≠ a) (called the node strength; NS) is given by Equation (11). Here, K represents the total number of electrodes K = 16.

For c₁ and |c₂|, repeated-measures analysis of variance (ANOVA) was performed to determine statistically significant differences between the younger and older groups. Age group was used as an inter-subject factor, and the 16 electrodes from Fp1 to T6 were used as intra-subject factors. The ANOVA results were represented by F values based on intra-group and inter-group variance comparisons. Greenhouse-Geisser adjustments were applied to the degrees of freedom. The α = 0.05 bilateral level was used; this was considered a statistically significant criterion for avoiding type I errors. A post-hoc t-test was subsequently used to evaluate the main effect between the younger and older age groups and effect of the interactions per electrode. Here, Benjamini-Hochberg false discovery rate (FDR) correction was applied to the t value for multiple comparisons of c₁ and |c₂| (q < 0.05) (16 p values: 16 electrodes).

In the ANOVA for sample entropy, age group was used as an inter-subject factor, and the 16 electrodes from Fp1 to T6 and a temporal scale were used as intra-subject factors. A post-hoc t-test was subsequently used to evaluate the main effect between the younger and older groups and effects of interaction for the electrodes and temporal scales. The α = 0.05 bilateral level was used. FDR correction was applied to the t scores for multiple comparisons (q < 0.05) (480 p values: 16 electrode × 30 scales).

In ANOVA for NS at each frequency band, age group was used as an inter-subject factor, and the 16 electrodes (from Fp1 to T6) were used as intra-subject factors. A post-hoc t-test was subsequently used to evaluate the main effect between the younger and older groups and effect of the interaction for the electrodes. The α = 0.05 bilateral level was used. FDR correction was applied to the t scores for multiple comparisons (q < 0.05) (80 p values: 16 electrodes × 5 frequency bands). The t-test was used for electrode-pair-wise group comparison of PLI between the younger and older groups. With a control for multiple comparisons, FDR correction was applied to the t scores (q < 0.05) (600 p values: 120 electrode pairs × 5 bands).

Older participants were classified using the receiver operating characteristics (ROC) curve. A logistic regression model based on sample entropy, c₁, |c₂|, and the NS of the PLI was used to identify older participants. Here, the logistic regression model outputs the “older participants” discrimination probability for each participant. The true/false positive rate at each threshold of discrimination probability from 0 to 1.0 for both groups was then measured. Principal component analysis was used as a pre-treatment for dimensionality reduction, and logistic regression based on c₁, |c₂|, sample entropy, and the NS of the PLI was implemented. The accuracy of discrimination was evaluated using the area under the ROC curve (AUC). We also used 5-fold cross-validation to prevent overfitting; AUC = 1.0 corresponds to perfect discrimination, and AUC = 0.5 corresponds to random discrimination. Here, the principal component analysis was conducted within cross validation (Shim et al., 2021) to avoid the inaccurate estimation of performance of discrimination. AUC values were averaged among 20 trials to choose tested and evaluated data set in 5-fold cross-validation and their standard deviations (SD) were also derived.

To evaluate the relationship between NS and complexity, we used Spearman's correlation R between NS and the complexity indexes (c₁, |c₂|). To control the multiple comparison, FDR correction was applied these R-scores (q < 0.05) (16 p values: 16 electrodes).

MF analysis was performed in the younger and older participants. Figure 2 shows the mean and standard deviation for each group with respect to D(h) and h. Owing to its wide distribution, this analysis is thought to reflect the multi-fractal characteristics (Sikdar et al., 2018) of the EEG signals for both groups. Table 1 shows the repeated-measures ANOVA results of the first (c₁) and second (|c₂|) cumulants within a singular spectrum. A strong main effects were observed for c₁ and |c₂|. The mean values of c₁ and |c₂| in the older and younger groups, as well as the results of a post-hoc t-test between the older and younger groups, are shown in Figure 3. The post-hoc t-test revealed that the value of c₁ was statistically significantly lower for the older participants at 13 electrodes (F3, Fz, F4, F7, F8, C3, C4, P3, Pz, P4, T6, O1, and O2) (q < 0.050). In addition, the results showed that the value of |c₂| was statistically significantly lower for the older participants at 14 electrodes (Fp1, Fp2, F3, Fz, F4, F7, F8, C3, C4, Pz, P4, T6, O1, and O2) (q < 0.050). The results of the MF analysis demonstrated that aging increases complexity (shown by less smoothness) and decreases multi-fractality. The time-series with large (small) multi-fractality exhibits intermittent and transient behavior with large (small) amplitude (Ihlen, 2012). Meanwhile, complexity reflects the degree of complexity for temporal behavior in entire the time-range, rather than intermittent behavior. Therefore, EEG signal in older subjects corresponds to homogeneous and highly complex temporal behaviors.

MSE analysis was performed in the younger and older participants. Table 2 shows the repeated-measures ANOVA results for the MSE analysis. We found that no main effect was observed, although there were interactions in the group × scale and the group × node × scale. The mean values of sample entropy in the older and younger groups, as well as the results of a post-hoc t-test between the older and younger groups, are shown in Figure 4. The results demonstrated a statistically significantly higher sample entropy for the older participants (q < 0.050) in the temporal-scale region of 1 to 5 (0.005 to 0.025 s at all electrodes). The results of the MSE analysis demonstrated that aging increases complexity on a fast temporal scale.

PLI analysis was performed on younger and older participants. Table 3 shows the ANOVA analysis results for the NS of the PLI for each band among younger and older participants. The results indicated that there was a main effect in the delta and gamma bands and that there was an interaction with respect to the group × node in the alpha, beta, and gamma bands. The post-hoc t-test results for the NS are shown in Figure 5C. Although no statistically significant differences satisfying with FDR criteria q < 0.05 were observed between the older and younger groups, relatively higher NS at delta and gamma band in the older group was observed. Regarding the PLI among pair-wise electrodes, the mean values of the PLI in the older and younger groups, as well as the results of t-tests between the older and younger groups, are shown in Figures 5A,B. No statistically significant differences satisfying with FDR criteria (q < 0.05) were observed between the older and younger groups.

To evaluate the relationship between complexity and functional connectivity, a correlation analysis was performed, using Spearman's correlation, between NS and the complexity indexes (c₁, |c₂|). Figure 6 shows the Spearman's correlations between NS of PLI and c₁ and between NS of PLI and |c₂| in both the younger and older groups. The correlation with c₁ did not meet the FDR correction criteria of q < 0.050; while there were positive correlations with NS at alpha and |c₂| at Fp1, Fp2, F3, F4, and Fz as well as NS at beta band and |c₂| at F3 in younger group. In Figure 7, the scatter plots at these electrodes were shown, significantly large positive correlations were observed. Therefore, large node strength might lead the intermittent and transient behavior reflecting |c₂| in local neural activity, instead of steady neural variability reflecting c₁.

To evaluate the classification ability for c₁, |c₂|, and the PLI, we evaluated the ROC. We observed a statistically significantly large sample entropy in the older group for time scales of 1 to 5, as shown in Figure 4. Therefore, the sample entropy was averaged in this temporal region for the purpose of clarification. Table 4 shows the results of the ROC in cases with the first-third principal components of c₁ and |c₂|, as well as sample entropy. In the results shown in Table 4, c₁ had the highest value (AUC = 0.86). Table 5 shows the results of the ROC in cases with the first-third principal components for the NS of the PLI in the delta, theta, alpha, beta, and gamma bands. In the results shown in Table 5, NS at alpha had the highest value (AUC = 0.84). Table 6 shows the results of the ROC in cases with the first-third principal components of the combination of c₁ and |c₂| and the combination of the NS at alpha, as well as c₁, |c₂| and sample entropy. Almost AUC values increased by combining the NS at alpha, c₁, |c₂|, and sample entropy in comparison with cases using a single index, as shown in Tables 4, 5. While, at the other bands, AUCs in the NS are significantly lower in comparison with the complexity indexes (c₁, c₂, and sample entropy) (see Tables 4, 5). Therefore, AUCs in the case with combinations of NS at the other band and the complexity index are inferior to AUCs in the case using a single complexity index.

To demonstrate that the determination area for older participants is determined by c₁, |c₂|, sample entropy, and the NS at alpha of the PLI, the determination area of the older participants was defined as the plane between the first principal components of c₁ and first principal components of |c₂|, and plane between the first principal components of c₁ and first principal components of the NS at alpha, plane between the first principal components of |c₂| and first principal components of the NS at alpha, and plane between the first principal components of sample entropy and first principal components of the NS at alpha (see Figure 8). Here, the other components, with the exception of the axis of the planes, were set to average among participants in both the younger and older study groups. The dependency on all of these factors in the decision region was confirmed.