Abstract

Understanding the temporal organization of brain activity requires methods that capture scale-free dynamics while accounting for the high-dimensional, spatially correlated nature of the electroencephalogram (EEG) data. We propose a novel framework that integrates nonlinear manifold learning (Isometric mapping) with detrended fluctuation analysis (DFA) to quantify long-range temporal correlations (LRTC) in the alpha-band of EEG signals. We applied this framework to two music related EEG datasets, as music is known to evoke different emotions and synchronize brain activity. The first dataset was obtained during live Indian classical music (ICM) listening that included two ragas, Yaman and Puriya Dhanashree. EEG was recorded from 13 healthy volunteers (24 channels, sampled at 500 Hz). The second dataset is the Music BCI dataset (006-2015), which includes Jazz and Synth-pop musical clips, with EEG collected from 11 subjects (64 channels, downsampled to 200 Hz). The EEG data from both datasets were preprocessed, band-limited to 8–13 Hz, and segmented into non-overlapping 2-s windows. Alpha-band power was extracted from each channel to form the feature matrix used for embedding. For the ICM dataset, Isometric mapping (Isomap) produced a low-dimensional representation (d = 3), which we analyzed using two approaches: (i) a norm-based approach and (ii) a mean-based approach. For comparison, an equivalent PCA-based pipeline (d = 5) was implemented. The Isomap mean-based DFA yielded consistent scaling exponents (α) in the range of 0.66–0.70, with higher goodness-of-fit (R2) and narrower bootstrap confidence intervals than the norm-based approach. PCA produced similar trends but required more dimensions. Paired t-tests showed that the Isomap mean-based approach detected music-related changes more sensitively than PCA (Yaman p = 0.02; Puriya Dhanashree p = 0.008). Comparable results were also observed for the second Music BCI dataset, where Isomap achieved a compact representation with d = 5, compared to d = 8 for PCA. In this dataset as well, the mean-based DFA yielded α values in the range of 0.62–0.65 and higher goodness-of-fit. Overall, the results suggest that combining nonlinear manifold embeddings with mean-based DFA provides a compact and robust framework for characterizing scale-free temporal structure in EEG data.

1 Introduction

Electroencephalography (EEG) records the brain's electrical activity from the scalp, yielding complex and nonlinear signals. Its superior temporal resolution is crucial for capturing rapid brain dynamics (Niedermeyer and da Silva, 2005). Numerous studies have shown that neural activity shows long-range temporal correlations (LRTC) consistent with fractal organization (Banerjee et al., 2016; Linkenkaer-Hansen et al., 2001). Assessing such correlations is important to understand how the brain maintains memory of past events and adapts to external stimuli (Hardstone et al., 2012). For quantifying LRTCs in nonstationary time series data, a well-known technique known as detrended fluctuation analysis (DFA) is used (Peng et al., 1995, 1994; Pavlov et al., 2020). DFA technique has been widely applied to EEG and magnetoencephalography (MEG) data to reveal scale-free dynamics during rest, task performance, and sleep (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012). These findings suggest that DFA can capture fundamental aspects of neural temporal organization across conditions, making it a suitable tool to investigate the effects of external stimuli such as music.

EEG signals are inherently nonlinear and exhibit strong spatial correlations across channels, resulting in high-dimensional yet redundant representations of brain activity (Stam, 2005; Sanei and Chambers, 2013; Nunez and Srinivasan, 2006). Such redundancy suggests that the underlying neural dynamics may lie on a lower-dimensional manifold embedded within the sensor space, as observed in studies of large-scale brain dynamics and electrophysiological recordings (Andrzejak et al., 2001; Mehrkanoon et al., 2014). Previous studies have shown that neural activity can often be characterized by low-dimensional representations, reflecting coordinated dynamics across brain regions (Mehrkanoon et al., 2014; Yu et al., 2022). To address this high-dimensional redundancy, dimensionality reduction techniques are used to extract informative components while preserving underlying dynamics. Linear dimensionality reduction techniques such as principal component analysis (PCA) are commonly used to obtain lower-dimensional representations (Jolliffe, 2011). The PCA captures the variance along orthogonal directions in Euclidean space, but fails to capture nonlinear dependencies. Whereas manifold learning methods are used to preserve the intrinsic geometry of the data by approximating geodesic distances on the underlying manifold. The Isometric mapping (Isomap) technique (Tenenbaum et al., 2000) is one of the nonlinear dimensionality reduction methods, which combines local nearest-neighbor graphs with multidimensional scaling (MDS) to recover global nonlinear structure. Applications of Isomap and related algorithms have demonstrated advantages in uncovering low-dimensional embeddings of complex biological and neural data (Saul and Roweis, 2003; Ashraf et al., 2023; Mitchell-Heggs et al., 2023; Anuragi et al., 2024). This motivates the use of nonlinear dimensionality reduction as a more robust representation of EEG dynamics compared to linear projections.

Beyond classical dimensionality reduction techniques, recent advances in deep learning have enabled powerful data-driven representation learning frameworks for EEG analysis. BCINetV1 (Aziz et al., 2025) employs a convolutional attention-based architecture that integrates temporal and spectral feature extraction to capture both local and global dependencies in non-stationary EEG signals. Similarly, large-scale models such as the Large Brain Model (LaBraM) (Jiang et al., 2024) utilize transformer-based architectures trained on multiple EEG datasets to learn generalized representations across diverse tasks. In addition, EEGPT (Wang et al., 2024) introduces a pretrained transformer model based on self-supervised learning, incorporating spatio-temporal representation alignment and reconstruction objectives to extract robust EEG features. While these approaches have demonstrated strong performance in EEG decoding and classification tasks, they primarily focus on learning hierarchical representations optimized for predictive performance. In contrast, the present work adopts a complementary perspective by focusing on interpretable and mathematically grounded analysis of EEG dynamics.

Despite the extensive use of DFA to quantify scale-free properties in EEG (Banerjee et al., 2016; Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012; Karkare et al., 2009; Grubov et al., 2025) and the growing application of dimensionality reduction techniques (Jolliffe, 2011; Tenenbaum et al., 2000; Saul and Roweis, 2003), relatively little work has integrated nonlinear manifold embeddings with fractal analysis. Most existing studies apply DFA directly to individual EEG channels, yielding interpretations restricted to those cortical regions (Banerjee et al., 2016; Pavlov et al., 2018) only. These approaches overlook broad synchrony across channels and fail to provide a global representation of brain dynamics. These limitations motivate the need for the development of a framework that combines nonlinear dimensionality reduction techniques with DFA to study long-range temporal correlations in EEG at a global scale. To the best of our knowledge, no prior literature works have explored this direction.

In contrast to traditional approaches where DFA is applied directly to EEG channels, we first reduce the dimension of the data using the Isomap technique. This approach assumes that the EEG data lie on an intrinsically curved, low-dimensional manifold rather than in the ambient sensor space. Isomap preserves this geodesic structure more effectively than linear methods. By applying DFA to the lower-dimensional embeddings obtained through Isomap, we characterize LRTCs, which more closely reflect the underlying brain dynamics. This integration of nonlinear dimensionality reduction with DFA provides a novel framework to uncover the hidden temporal structures in EEG activity. The main contributions of this work are as follows: Firstly, we propose a framework that integrates Isomap with detrended fluctuation analysis (DFA) for the study of EEG dynamics. Secondly, within this framework, we opted for two approaches to perform DFA on low-dimensional embeddings: a norm-based approach and a mean-based approach. Thirdly, we compare the proposed framework with the traditional linear dimensionality reduction technique using PCA.

2 Dataset and subjects2.1 Subjects

For this study, 13 healthy individuals were recruited, including 11 males and two females, with a mean age of 30.54 ± 9.42 years. The demographic details of the participants are provided in Table 1. None of the participants had a history of neurological disorders and hearing loss. Before data collection, a clear demonstration of the experimental procedure was given to each subject. EEG data was recorded using a 24-channel cap, with electrodes placed according to the standard international 10–20 system (Jasper, 1958). The EEG data were acquired with the CameraEEG application using a mobile phone (Hazarika et al., 2023). To ensure signal quality, electrodes were cleaned with alcohol, conductive gel was applied, and impedances were maintained between 10 and 20 kΩ. Permission was obtained from the institute's ethics committee at IIT Guwahati for collecting data from human participants.

SubjectsGenderAge (years)Prior music knowledgeSubject 1M21NoSubject 2M45NoSubject 3M26NoSubject 4F25NoSubject 5M31NoSubject 6M28NoSubject 7M28YesSubject 8M41YesSubject 9M27NoSubject 10M25NoSubject 11F52YesSubject 12M25NoSubject 13M23Yes

Demographic data of the volunteers participated in the study.

2.2 Dataset used

We used an EEG dataset from a music paradigm employed in the studies (Panda et al., 2024, 2023). The EEG data were obtained while subjects listened to a live vocal musical paradigm. It consisted of two ragas: Yaman and Puriya Dhanashree (PD). Data were collected on different days depending on the availability of participants. The musical paradigm is shown in Figure 1. It begins with an initiation phase, which lasts 30–45 s, and then enters Raga Yaman. This is followed by a relaxation period and then a transition phase. After the transition phase, raga Puriya Dhanashree (PD) is played, followed by another relaxation period. No music was played during either of the relaxation periods, although the EEG data were recorded. Each raga (Yaman and PD) lasted 3 min, and both relaxation conditions were 2 min long. The transition phase comprised 1 min of Raga Yaman followed by 1 min of Raga PD. Since it was a live musical paradigm, the timestamps of music and relaxation periods varied across subjects.

Live musical paradigm used for data collection for this study. Consists of two ragas: Yaman and Puriya Dhanashree, and relaxation conditions.

2.3 Preprocessing

EEG data were recorded while participants listened to the musical paradigm with their eyes closed. During the relaxation periods, participants were allowed to open their eyes and observe the surroundings. The data were sampled at a frequency of 500 Hz. Faulty channels were excluded from the study, and then artifact correction was performed using the multichannel wiener filter (MWF) toolbox in MATLAB (Somers et al., 2018). It removes eye blinks, muscle, and movement artifacts. Then, the EEG data were normalized using z-score normalization. Finally, the EEG signals were band-pass filtered in the alpha band (8–13 Hz). All subsequent analyses in this study were conducted on the alpha band (Panda et al., 2024).

3 Methodology

This section outlines the methodology adopted for applying detrended fluctuation analysis (DFA) on the low-dimensional manifold obtained using the isometric mapping (Isomap) technique. The pipeline consists of three main steps: (i). Extraction of alpha band (8–13 Hz) power, (ii). Dimensionality reduction of EEG data using Isomap, and (iii). Application of the DFA algorithm to the low-dimensional embedding.

3.1 Alpha-band (8–13 Hz) power

EEG data were acquired using a 24-channel cap at a sampling frequency (fs) of 500 Hz. The recordings are represented by a matrix X∈ℝN×T, where N is the number of channels (N = 24) and T the number of time points. Then, the EEG signal was segmented into n non-overlapping 2-s windows. The number of windows n was computed using the ‘floor' function in MATLAB. This ensures that only complete windows are considered. Any remaining samples at the end of the recording that do not form a full window were discarded and not included in the analysis.

In Equation 1, the window length (L = 2·fs) of 2 s was chosen to provide sufficient frequency resolution (0.5 Hz) for capturing the alpha band of the EEG signal, consistent with prior EEG studies (Klimesch, 1999; Sanei and Chambers, 2013). Although EEG signals may not be strictly stationary within these windows, detrended fluctuation analysis (DFA) is specifically designed to handle nonstationary time series. The detrending step in DFA removes local trends within each segment, thereby reducing the influence of nonstationarities on the estimation of the scaling exponents (Linkenkaer-Hansen et al., 2001; Peng et al., 1994; Kantelhardt et al., 2001). For each 2-s window, the alpha-band power was computed for every channel using MATLAB's function ‘bandpower'. This function estimates the power spectral density (PSD) using Welch's method (Fast Fourier Transform-based) and then integrates the power within the alpha-band (8–13 Hz). These values were stored in a matrix , where n is the number of windows and N the number of channels (N = 24). For subjects with one or two faulty channels removed, N = 24 was reduced to 23 or 22, respectively.

3.2 Dimensionality reduction of EEG data using Isomap

To reduce the dimensionality of the EEG data matrix Xα, the nonlinear manifold learning technique: Isomap, was employed. Unlike linear methods, Isomap preserves the intrinsic geometry of the EEG data by maintaining geodesic distances between points on the underlying manifold. This algorithm consists of three major steps (Tenenbaum et al., 2000):

3.2.1 Construction of neighborhood graph

To preserve the local geometry of the EEG data, the k-nearest neighbors (k-NN) approach is employed. In this step, a graph is constructed by connecting each segment (time window) to its k nearest neighbors using the Euclidean distance. The pair wise Euclidean distance is calculated using Xα. The pairwise Euclidean distance between time windows i and j is:

where Xi and Xj in Equation 2 denotes the vectors corresponding to time windows i and j, respectively. The distance matrix D∈ℝn×n contains all pairwise distances. For each point i, its k-nearest neighbors Nk(i) are determined by selecting the k smallest entries in the distance matrix D. Using this, an adjacency matrix W is constructed using Equation 3.

The resulting adjacency matrix W∈ℝn×n encodes neighborhood relations among the n time windows: Wij < ∞ when window j is among the k-nearest neighbors of i (or vice versa), and Wij = ∞ otherwise. The connectivity of the resulting graph is then verified, and the optimal value of k is selected such that the graph becomes fully connected. The procedure for determining this optimal k is outlined in Algorithm 1.

Optimal value selection.

In Algorithm 1, the initial value is chosen as as a heuristic initialization (Zhang, 2016). The symbols V and E denote the set of vertices (time windows) and edges (connections) of the graph, respectively. The neighborhood size k is selected adaptively rather than being fixed. Starting from k0, the value of k was incrementally increased until the k-nearest neighbor graph becomes fully connected. Among all such values, the smallest k ensuring connectivity is selected.

3.2.2 Estimation of geodesic distances

The goal of this subsection is to find the true geodesic distances by computing the shortest paths on the neighborhood graph (W). The geodesic distance matrix (G) is computed using Dijkstra's algorithm (Dijkstra, 1959, 2022). This matrix G∈ℝn×n represents the geodesic distances between all pairs of time windows.

3.2.3 Construction of low-dimensional embedding

In this step, the aim is to obtain a low-dimensional representation of the EEG data, Y∈ℝn×d, where d is the required lower dimensions while preserving the geodesic distances. Classical multidimensional scaling (MDS) is performed on the geodesic distance matrix (G) to achieve this. Compute the kernel matrix (B) using the Equation 4.

where , and H is the centering matrix, given by

The identity matrix (I) in Equation 5 is of size n×n, and 11⊤ denotes the matrix of all-ones, and of size n×n. By performing the eigen decomposition of the kernel matrix B, we obtain

where Λ and V in Equation 6 denote the eigenvalues and eigenvectors of B, respectively. Select the top d eigenvalues and eigenvectors, and form the embedding matrix Y as,

From Equation 7, the original EEG data is reduced to a lower-dimensional embedding matrix Y∈ℝn×d. The complete algorithm of Isomap is presented in Algorithm 2.

Isomap dimensionality reduction.

3.3 Application of DFA algorithm

From Equation 7, the low-dimensional embedding matrix Y∈ℝn×d represents the global structure of the EEG data in a reduced d-dimensional space. Two approaches can be considered for the DFA analysis: (i) applying DFA to a one-dimensional signal obtained by taking the Euclidean norm across the d dimensions, (ii) applying DFA separately to each dimension of Y and averaging the resulting scaling exponents (α).

The mathematical formulation below is presented for the norm-based approach, while the detailed formulation of the mean-based approach is provided in the Appendix.

For applying the proposed DFA algorithm, the low-dimensional embedding matrix Y∈ℝn×d was converted into a one-dimensional representation of size n×1 by computing the Euclidean norm across the d dimensions.

Y∈ℝn×d in Equation 8 is the low-dimensional embedding, and is the resulting one-dimensional time series. This is then used as the input for DFA.

The DFA technique has three major steps (Peng et al., 1994). The first step is, mean removal and cumulative sum. For the given signal Ynorm(i), i = 1, 2, …, n, the mean is defined as

The cumulative sum after removal of the mean is computed using Equation 9 as,

The second step of the DFA algorithm is the local detrending, in which, the integrated series from Equation 10Z(k) is divided into p non-overlapping segments of equal length. It is varied within the range of,

Within each segment, a least squares polynomial fit Zp(k) is computed to represent the local trend. The root-mean-square fluctuation at scale p of the integrated series is given by,

Z(k)−Zp(k) in Equation 12 is known as detrending. The range in Equation 11, [4, n/4] (Peng et al., 1994) ensures that the analysis includes multiple scales while maintaining statistical reliability, as both very small and very large box sizes are avoided.

The final step of DFA is the scaling behavior. Now, the relationship between the detrended series and the segment lengths given as,

The α in the Equation 13 is known as the scaling exponent and is defined as the slope of the double logarithmic plot log F(p) vs. log(p). The complete DFA on the low-dimensional embedding is presented in Algorithm 3. The scaling exponent (α) provides the long-range temporal correlation (LRTC) of the EEG data. When α = 0.5, it indicates the EEG signal is completely uncorrelated and represents white noise. When long-range temporal correlations are analyzed in EEG signals, a power law pattern typically produces scaling exponents between 0.5 and 1. Values closer to 1 indicate stronger persistence, indicating the correlations decay more slowly over time. If the exponent exceeds 1, the relationship no longer follows a power law, and LRTC is not preserved (Linkenkaer-Hansen et al., 2001). The complete bird's-eye view of the proposed methodology is shown in Figure 2.

Detrended fluctuation analysis (DFA).

Graphical abstract of the proposed framework showing the preprocessing, Isomap embedding, and DFA analysis steps.

4 Results4.1 Selection of lower-dimensional embedding

Finding the appropriate number of lower dimensions after constructing the embedding matrix (Equation 7) is crucial for the proposed DFA analysis. To determine the optimal embedding dimension, the elbow method was employed on the eigenvalue spectrum of the kernel matrix B (Equation 6). This method identifies the point where the eigenvalue decay curve exhibits a clear inflection, suggesting, the additional dimensions contribute only marginally to the variance.

Figure 3 illustrates the eigenvalue spectra for both ragas in subject 12. In both ragas, the curve exhibits a clear elbow after the third dimension, beyond which the eigenvalues level off. Based on this observation, we have selected three (3) dimensions for the reduced embedding. Figures 3a, b demonstrate the eigenvalue spectrum behavior for both the ragas (Yaman and Puriya Dhanashree, respectively). The similarity in the eigenvalue spectral decay in both ragas supports the robustness of choosing a three-dimensional embedding. This ensures the essential geometric structure of the EEG data is preserved while minimizing redundancy and noise. While Figure 3 shows the eigenvalue spectrum for one subject, consistent spectral patterns were observed in all subjects and both ragas (refer to Supplementary File 1). This consistency supports the adoption of d = 3 as the embedding dimension throughout.

Eigenvalue spectra of the kernel matrix (B) for Subject 12. The elbow after the third eigenvalue indicates an optimal embedding dimension of d = 3. (a) Yaman. (b) Puriya Dhanashree.

To provide a more quantitative comparison, we additionally evaluated neighborhood preservation using the trustworthiness metric (Venna and Kaski, 2001) and variance explained (Tables 2, 3). From Table 2, it is observed that PCA shows higher trustworthiness values (≈0.98–0.99), which reflects the preservation of local Euclidean structure. In contrast, Isomap achieves consistently high trustworthiness (≈0.94–0.96), indicating robust neighborhood preservation while additionally capturing nonlinear geometric structure. The trustworthiness metric was computed using neighborhood sizes Kt = 5 and Kt = 7, where Kt denotes the number of nearest neighbors considered for the local structure preservation. These values were selected based on prior studies (Venna and Kaski, 2001; Van Der Maaten et al., 2009), where small neighborhood sizes (5–15) are commonly used to assess local neighborhood preservation. Furthermore, the three-dimensional (d = 3) embedding captures a substantial portion of the variance, with mean explained variance of 80.02%±8.12% for Raga Yaman and 82.30%±8.96% for Raga Puriya Dhanashree. For nonlinear methods such as Isomap, variance explained alone is not sufficient to characterize embedding quality, as the method is designed to preserve intrinsic geometric structure rather than maximize variance.

MethodKtYamanPuriya DhanashreeDurinAfterDuringAfterIsomap50.957 ± 0.0150.937 ± 0.0340.956 ± 0.0140.935 ± 0.03270.961 ± 0.0140.944 ± 0.0300.960 ± 0.0140.936 ± 0.034PCA50.993 ± 0.0050.984 ± 0.0160.993 ± 0.0040.984 ± 0.01270.994 ± 0.0040.986 ± 0.0140.994 ± 0.0050.986 ± 0.010

Trustworthiness (Mean ± Std) for both the methods PCA and Isomap during and after music listening for Raga Yaman and Raga Puriya Dhanashree.

SubjectYaman (%)Puriya Dhanashree (%)180.1579.40284.8697.47382.7498.84477.0276.70574.5588.21685.3986.65789.5986.79877.3380.83973.2274.061060.4772.111186.8475.011289.6285.291378.5472.54

Cumulative explained variance (%) captured by the first three dimensions (d = 3) for all subjects across Raga Yaman and Raga Puriya Dhanashree.

4.2 Temporal scaling exponent (α)

The temporal scaling exponent (α) characterizes the strength of long-range temporal correlations (LRTC) in the EEG dynamics.

The deviations of α from the white noise baseline quantify the degree to which the EEG exhibits temporal self-similarity. Figures 4, 5 illustrate the scaling behavior by showing the log–log dependence of the fluctuation function F(p) on the scale p for both the ragas Yaman and Puriya Dhanashree, respectively. For each raga, the results are presented for both the music and relaxation conditions. The linear trends demonstrate the presence of scale-free dynamics, while the slopes of the fitted lines yield the scaling exponents (α). The differences between the music and relaxation conditions are reflected in distinct α values, indicating that musical stimulation modulates the persistence of temporal correlations in EEG activity.