Complexity in Speech and Music Listening via Neural
Manifold Flows

Claudio Runfola¹*, Matteo Neri¹,², Daniele Schön¹, Benjamin Morillon¹, Giovanni Rabuffo¹, Pierpaolo Sorrentino¹* and Viktor Jirsa¹

1. Aix Marseille Université, INSERM, INS, Institut de Neurosciences des Systèmes, Marseille, France

2. Aix Marseille Université, CNRS, INT, Institut de Neurosciences de la Timone, Marseille, France

senior authors

*corresponding authors:

claudio.runfola@univ-amu.fr, pierpaolo.sorrentino@univ-amu.fr

Abstract

Understanding the complex neural mechanisms underlying speech and music perception remains a multifaceted challenge. In this study, we investigated neural dynamics using human intracranial recordings. Employing a novel approach based on low-dimensional reduction techniques, the Manifold Density Flow (MDF), we quantified the complexity of brain dynamics during naturalistic speech and music listening and during resting state. Our results reveal higher complexity in patterns of interdependence between different brain regions during speech and music listening compared to rest, suggesting that the cognitive demands of speech and music listening drive the brain dynamics towards states not observed during rest. Moreover, speech listening has more complexity than music, highlighting the nuanced differences in cognitive demands between these two auditory domains. Additionally, we validated the efficacy of the MDF method through experimentation on a toy model and compared its effectiveness in capturing the complexity of brain dynamics induced by cognitive tasks with another established technique in the literature. Overall, our findings provide a new method to quantify the complexity of brain activity by studying its temporal evolution on a low dimensional manifold, suggesting insights that are invisible to traditional methodologies in the contexts of speech and music perception.
# Introduction

Music and speech processing are among the most complex, characteristically human, cognitive functions. They rely on the modulation of acoustic parameters to yield abstract representations of sound sequences, have a large range of time constants, and significant information. In both, complexity is built up following a hierarchical organization by which a finite set of elements (phonemes or notes) are combined into more complex entities (words and melodies), which in turn are organized to generate more complex ensembles (such as sentences or refrains), arranged across several other levels to produce natural speech and music (1,2). How the human brain attributes meaning and qualities to sounds remains an open and intricate question (3–6).

A long tradition in neuroscience supports the hypothesis that the interactions between neural units at different spatiotemporal scales are fundamental for sustaining complex cognitive functions and interpreting such stimuli (7,8). Complexity arising from finely-tuned interactions among neuronal populations may be mirrored in low-dimensional patterns of interdependence among neural signals. To study these, a great deal of success was obtained by methodologies and theoretical approaches that allow a low-dimensional representation of multivariate brain recordings (9). A growing body of evidence demonstrates that the dynamics of the interactions plays a crucial role in sustaining cognitive tasks (10–12), including speech (13–18) and music (19) perception.

However, quantifying and comparing the complexity of the temporal evolution of low-dimensional brain dynamics underlying speech or music tasks (beyond focusing on topological properties and task-related segmentation of the latent space) remains challenging. Furthermore, most previous works employed low-dimensional descriptions of brain activity to investigate simple tasks which allow the isolation of cognitive features of interest in controlled scenarios, but provide limited ecological validity (13–17). Moreover, while short and simple stimuli in the case of speech might still induce linguistic processing, they are often not cognitively relevant to music listening (4,5,20). Another relevant point concerns brain recording techniques. As an example, functional magnetic resonance (fMRI) fails to capture the fast dynamics of music and speech due to its low temporal resolution (17).

To address these points, we leveraged a data set in which the brain activity of 19 epileptic patients were recorded using stereotactic electroencephalography (sEEG) during naturalistic speech and music listening (Fig. 1). The ecologically-relevant
stimuli provide a genuine perspective on speech and music listening complexity in a
naturalistic scenario. Furthermore, the sEEG signal provides anatomically precise
information about the functionally selective engagement of neuronal populations at
the millimeter scale and the temporal dynamics of their engagement at the
millisecond scale, which is necessary for depicting accurately the neurophysiological
underpinning of specific cognitive processes (21,22). We designed a novel
methodology aimed at mapping the dynamical changes underlying the system's
evolution in a neural manifold across diverse tasks, enabling quantitative and
comparative analyses of complexity. This method accounts for probabilistic density
evolution within low-dimensional spaces defined by latent variables, offering
complementarity to existing techniques for assessing the topological properties and
task-related segmentation of the manifolds onto which high-dimensional neuronal
activity collapses (23–31).

We hypothesize that speech and music listening display more complex large-scale
dynamics, as compared to rest. The hypothesis is based on the more complex
cognitive demands that music and speech listening impose on the brain. To further
clarify the interpretation of our findings, we first investigate a simple in-silico model
with known dynamics.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Brain activity was recorded using stereotactic electroencephalography (sEEG). Participants passively listened to a 10-minute storytelling segment and a 10-minute music piece, both embedded within three resting-state sessions (each lasting 3.3 minutes). Latent system evolution during the different tasks is captured by a low-dimensional unfolding of the neural space. Comparing dynamical changes across neural manifolds enables a comparative analysis of complexity.</figcaption></figure>
Results

For simulated or empirical data, we performed an independent component analysis
(ICA) on the time series and we projected their time evolution on two components
(Fig. 2A). Then, the density maps of the system's trajectory in different time windows
(see Fig. 2B and Methods section) were correlated between one another to obtain
the Manifold Density Flow (MDF) matrix (Fig.2C). In this way the element MDF(t₁, t₂)
corresponds to the correlation between the density maps of the system trajectory in
the space of the two components of the ICA, during the time-window w₁ (starting at
instant of time t₁) and the time-window w₂. In other words, MDF(t₁, t₂) measures the
extent to which the exploration of the latent space in w₁ and w₂ is similar. Considering
all the time windows allows us to quantify the complexity of the system's evolution, in
terms of the thoroughness of the exploration of the manifold onto which it is
embedded. Estimating the extent to which the manifold is (un)stable across time is a
direct measure of its complexity. Intuitively, more complex dynamics, characterized
by a greater number of attractive states, evolve more smoothly onto a
low-dimensional space, making large jumps less frequent; this behavior is reflected
in lower variance of the MDF variance. This information about the evolution of the
dynamics could not be unambiguously captured by the mean of the MDF, since
different types of dynamics might produce the same mean MDF. Note that, given the
symmetric nature of the MDF matrix, i.e. MDF(t₁,t₂)=MDF(t₂, t₁), we optimize our
calculations by focusing only the upper triangular part of the matrix.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Schematic picture of the construction of the manifold density flow (MDF) matrix. (A) Through dimensionality reduction techniques, signals extracted from high-dimensional datasets are projected onto a two-dimensional space in which the underlying dynamics of the system takes place. (B) Time-series of the two obtained components are segmented into overlapping time windows using a sliding window approach. For each time window, a probability density map is estimated from the distribution of the components in the previously obtained two-dimensional space. Despite the complex manifold architecture spanned by the independent components, simpler patterns emerge in the probability distribution flow, corresponding to the transitions between different system states. (C) Computing the Pearson's correlation between each pair of these probability densities yields the MDF matrix. Each entry *ij* quantifies the degree of similarity between the probability density map at time-window *w*<sub>*i*</sub> and that at time-window *w*<sub>*j*</sub>. The MDF captures the similarity in exploration of the latent space between different segments and allows for the extraction of quantifiable metrics of the system's inner dynamics complexity.</figcaption></figure>

**Simulations of known dynamics.** We first utilized a theoretical model with known dynamical complexity to test this interpretation of our method. We simulated a dynamical system composed of *N*=2,3,4 independent variables (Fig. 3A). The time evolution of each variable was produced by a cubic differential equation (see Methods). Noise allows for intermittent transition between a down- and an up-states, corresponding to two equilibrium points, with more noise resulting in more frequent transitions. Note that the number of attractive states scales with the number of variables, as observed from the projection of the system into the two ICAs
components (Fig. 3B). A total of 9 scenarios were simulated by manipulating the number of variables and the amount of noise, the latter directly affecting the switching rate between the different states of the system. With the simulated dynamics, we computed the MDF matrix using overlapping windows with a length of 800 time points and a slide of 10 time points. We observed significantly lower standard deviation in the MDF either when increasing the amount of noise, i.e. the switching rate, or the number of independent variables, i.e. the number of states explored by the system (Fig. 3C).

**During speech and music listening brain activity explores richer dynamical patterns.** Our results on synthetic data suggest that a decreasing standard deviation of the MDF could be due either to an increase in the number of states explored by brain dynamics or to a higher noise-induced switching rate. To assess whether there is a significant difference in the switching rate, we binarized the data using different thresholds (2 to 3 standard deviations; step=0.1). Then, we counted the number of 'switches' for each brain region, i.e. how often the region crossed the threshold, and averaged across regions for each participant and condition. For all the thresholds tested, we did not find any significant difference between speech, music, and rest (Wilcoxon test: all P>0.05) (Fig. 3D, left). Following this, we computed the MDF matrix during the three conditions using overlapping windows with lengths varying from 1 to 30 seconds and a slide of 100 ms. The standard deviation of the MDF during both speech and music listening is significantly lower than during rest (Cluster-based permutation test (CBPT): P=0.005 and P=0.026, respectively; Fig. 3D, center and right), with the standard deviation of the MDF during music listening being significantly lower than during speech (P=0.033). This supports the hypothesis that brain activity explores more states during speech than music listening, and fewer during rest. This interpretation suggests that complex cognitive functions need to be sustained by correspondingly more intricate brain dynamics. Furthermore, to assess the efficacy of our methodology, we conducted a comparative analysis between our results and those derived from the use of the standard deviation of the dynamic Functional Connectivity (dFC) matrix (32) as a metric to quantify and compare complexity across the three distinct conditions. We show how the use of this metric fails to tell the three conditions apart (CBPT: P>0.07)(Supplementary Material). Notably, our findings are confirmed when using three independent components instead of two, as detailed in the Supplementary Material (S2).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top panels: **Simulations.** From left to right the simulated systems with an increasing number of state variables (2,3 and 4). (A) Separated time series, represented in a carpet plot, of the bistable variables. (B) Representation of the time series in the space of the two components of the ICA, alongside with the upper triangular part of the MDF. (C) Distributions, extracted from 100 samples (100 simulations performed per each condition), of the standard deviations σ of the MDF matrices. The colors indicate varying levels of noise, which correlates with the switching-rate, ranging from low (green) to high (purple). Variance in the MDF matrix decreases either when increasing in the number of variables or the amplitude of noise, mirroring the increase in the number of states arising from the superposition of bistable regimes and the switching between them. Bottom panels: **Real data.** (D) On the left are the distributions of the switching rates, computed averaging the number of switches after having binarized the data (threshold=2.5 standard deviations); no significant differences are observed between speech, music, and rest conditions (Wilcoxon test: P>0.05). At the center, the means of the distributions of the standard deviation σ of the MDF matrices plotted against different numbers of time-window length (in seconds) employed in the computation of the MDF. Speech and music listening display higher complexity than rest (CBPT: P=0.005 and P=0.026, respectively), while music listening maintains simpler dynamics compared to speech (P=0.033). This behavior is more clearly illustrated on the right, where we present the distributions for a 15-second window length used to compute the MDF, along with the upper triangular sections of the MDFs for one participant across the three different tasks.</figcaption></figure>
**Brain dynamics is more synchronous during speech listening.** We conducted an inter-subject correlation analysis to evaluate the consistency of brain activity patterns across participants for the three different tasks. Specifically, we assessed Pearson's correlation between the MDF matrices computed for each subject. We observed significantly higher inter-subject correlations during speech listening compared to both music listening (CBPT: P=0.0001) and resting state (CBPT: P=0.0001) (Fig. 4).

<figure><img src="image_5.png" /><figcaption>**Figure 4. Synchrony in brain dynamical patterns.** On the left, the mean correlation matrices ($\rho$) between the MDF matrices are plotted against various time-window lengths (in seconds) used in the MDF computation. Speech listening exhibits higher synchrony in brain activity compared to music listening and rest (CBPT: P=0.0001 for both comparisons). On the right, the inter-subject correlation matrices for speech listening, music listening, and rest are shown, calculated using a 15-second time window for the MDF computation; a more uniform pattern of synchronization could be observed during speech listening.</figcaption></figure>

Discussion

In this paper, we set out to test the hypothesis that passive speech and music listening demands a more complex interplay among brain regions, as compared to resting state. Such interactions are mirrored in intricate patterns of interdependencies among recorded brain signals. To quantify these, we analyzed a data set collected using stereotactic intracranial electroencephalography to measure brain activity of 19 epileptic patients during naturalistic speech and music listening and during resting state. According to our hypothesis, we found a growing richness of the dynamical repertoire from rest to music to speech, as quantified by the corresponding diminishing standard deviation of the Manifold Density Flow (MDF), which expresses the ability to occupy particular dynamical states and transition
between them. This is in accordance with previous literature showing that the
presence of tasks pushes the brain dynamics towards configurations that would be
otherwise inaccessible during resting state (33).

Our results can be framed within the literature that uses dimensionality reduction
techniques to distinguish the dynamics that are induced by different natural stimuli
(10,33). However, somewhat differently than with fMRI, dimensionality reduction is
notoriously challenging when dealing with neurophysiological data. Classical
manifold approaches focused on the classification of the states constituting the latent
space, while disregarding the information about the flow among such states (13–17).
To overcome these challenges, our paper proposed an approach to characterize the
evolution of the density maps obtained from dimensionality reduction over time,
successfully capturing the different dynamical complexity in response to different
naturalistic stimuli. Our approach is influenced by choices such as the time-window
duration. We show that our results are generally stable over different windows'
durations and we highlight that these choices allow us to characterize and
investigate different temporal scales.

Of note, we measure the complexity of the dynamics in terms of the standard
deviation of the MDF. This metric captures changes in the system's dynamic
patterns, reflected by shifts in trajectories and transitions between attractors in the
phase space. This way, it might serve as a quantitative and comparative measure of
the speed of transitions between brain dynamical states, accounting for the diversity
of states (or attractors in the neural manifold) and the likelihood of transitioning from
one to another. In our toy model, a state is defined as the combination of active and
non-active units and we show in a synthetic model that a higher number of possible
visited states corresponds to a lower standard deviation of the MDF. One challenge
is that higher noise might artificially increase the complexity of state space
exploration, simply by increasing the probability of any region changing state and so
the switching rate. Hence, we show that, keeping the noise level constant, increasing
the number of “true” dimensions of the system, i.e. its complexity, unambiguously
results in lowering the standard deviation of the MDF. In analogy, we have checked
that the switching rate of the channels in the real data (that is, the average number of
switches across regions) does not change among conditions, supporting a qualitative
change in the dynamics. The comparison of the MDF to the dynamic Functional
Connectivity (dFC) finally shows that the MDF metric emerges as a more suitable
tool for analyzing the complex dynamics stemming from regions with higher
integration levels. Moreover, our results do not directly depend on the number of ICA
components chosen, as demonstrated by the same trends observed using three
independent components instead of two (see Supplementary Material (S2)).

We performed an inter-subject correlation analysis to evaluate how synchronously
brain dynamics is expressed among different individuals during the three cognitive
tasks. We observed that brain dynamics during speech listening is more
synchronous than during music listening or rest, indicating a common
spatio-temporal structure imposed by the external stimulus. This finding is consistent
across various time-window lengths used to compute the MDF matrix and further
underscore the utility of the MDF in capturing common neural responses among participants, particularly in the presence of a common external stimulus.

In conclusion, in this paper, we used intracranial data recorded during naturalistic speech and music listening and advanced dimensionality reduction techniques to study how the evolution of the interactions among brain regions is modulated by the cognitive conditions. We demonstrated the growing complexity of the dynamics from rest to music to speech. These may reflect the complexity of the neural mechanisms sustaining information integration and segregation during these different cognitive conditions. Future analysis may include a more detailed characterization of the regions involved in these processes.

# Materials and methods

**The model.** The system in Results section (a) is described by the node equation

$$x_i = (1 - x_i^2)x_i \quad (1)$$

in which the nodes are allowed to have intrinsic bistable dynamics. Integration is performed through the stochastic differential equation (SDE) scheme

$$\psi(t) = f(\psi(t))dt + g(\psi(t)) \otimes dW(t),$$

in which $\psi$, $f$, $g$ and $W$ are vectors of same dimension and $\otimes$ indicate the component-wise product. In particular, $\psi$ is the state vector $\psi = [\psi_1, \dots, \psi_N]$, with $\psi_i = x_i$, $f$ is the vector whose components are derived from the function on the r.h.s. of (1), i.e. $f_i = x_i$, $g = [g_1, \dots, g_N]$, and $dW(t)$ is a differential of a Wiener process with Gaussian increment. All the nodes are set to have the same level of Gaussian noise $g_i = \sqrt{2\sigma}$, where $\sigma$ tunes the noise amplitude. The system is simulated using just-in-time compilation for stochastic differential equations (JITCSDE) (34), a Python implementation of the adaptive integration method proposed by Rackauckas and Nie (35), with a sampling rate $dt=0.01$ over the interval $[0,2000]$, resulting in time series comprised of 200000 time points.
**The Manifold Density Flow matrix.** We make use of Independent Component Analysis (ICA) to embed the continuous time series data into a two-dimensional space. Subsequently, we segment the resulting time series of the two independent components into smaller time windows of length *L*, employing a sliding window approach where slides occur every *s* time points. Within each time window *T*ᵢ, we generate a smoothed 2D probability density estimation, PDE(*T*ᵢ), using the values of first and second independent components extracted from the time points spanning the interval starting at *T*ᵢ and ending at *T*ᵢ+*L*. This involves constructing a 2D histogram of the data using a specified number of bins along both dimensions, capturing the distribution of the data points within the time window interval. For each time window, the probability density is estimated using the boundaries defined by the space spanned by the two independent components across all time points. To achieve a smoother estimate of the probability density, we apply a Gaussian filter with a specified standard deviation to the generated histogram. This filtering process reduces noise in the density estimate, resulting in a smoother representation of the underlying distribution. In our analysis, we employ 128 bins for the 2D histogram and a standard deviation s.d.=8 for the Gaussian filter. Subsequently, we calculate the Pearson's correlation between each pair of these probability densities (heatmaps) to characterize dynamic changes in the density of the manifold generated by the independent components over time. This correlation analysis yields to a matrix, the Manifold Density Flow (MDF), which encapsulates the evolving relationships between the probability density distributions on the low-dimensional space across different time windows. The procedure we propose can be extended, in principle, to an arbitrary number of components *N*. In such cases, the correlation should be computed between ND density maps to generate the MDF matrix. In Supplementary Material (S2), we demonstrate how this can be done using three independent components as an example.

**Inter-subject correlation analysis.** MDF matrices were computed for each participant under the three different cognitive conditions and across various time-window lengths. To facilitate comparison, the MDF matrices were flattened into vectors. Pearson correlation coefficients were then calculated between these flattened MDF vectors for different participants, resulting in a correlation matrix that quantifies the synchrony in brain response to the cognitive stimulus across subjects.

**Participants.** The study involved 19 participants, comprising 10 females, with a mean age of 30 years (range: 8-54 years), all diagnosed with pharmacoresistant epilepsy. All participants were native French speakers. Before undergoing stereotactic intracranial EEG (sEEG) recordings, neuropsychological assessments were conducted to ensure intact language functions and normal hearing, as per the criteria. Notably, none of the participants had their auditory areas identified as part of the epileptogenic zone, confirmed by an experienced epileptologist's assessment. The recordings were conducted at Hôpital de La Timone in Marseille, France. Prior
to the experimental session, participants provided informed consent, and the
experimental protocol received approval from the Institutional Review Board of the
French Institute of Health (IRB00003888).

**Data acquisition.** The sEEG signal was captured utilizing depth electrode shafts with a diameter of 0.8 mm, containing 10 to 15 electrode contacts manufactured by Dixi Medical or Alcis from Besançon, France. Each electrode contact was 2 mm long and spaced 1.5 mm apart from each other along the shaft. Electrode implantation locations were determined solely based on clinical considerations. Participants were included if their electrode map partially covered Heschl's gyrus, either on the left or right side. The cohort comprised 13 unilateral implantations (10 left, 3 right) and 6 bilateral implantations, totaling 271 electrodes and 3371 contacts. During recordings, patients were situated either in an insulated Faraday cage or in a bedroom setting. In the Faraday cage, participants were comfortably seated in a chair within a sound-attenuated room. Data acquisition was performed using a 256-channels amplifier manufactured by Brain Products, with a sampling rate of 1 kHz and a high-pass filter set at 0.016 Hz. Alternatively, in the bedroom setup, data were collected using a 256-channels Natus amplifier from the Deltamed system, sampled at 512 Hz, and high-pass filtered at 0.16 Hz.

**Experimental design.** Participants in the study passively listened to a 10-minute storytelling segment (576.7 seconds, "La sorcière de la rue Mouffetard," Gripari P. 2004) and a 10-minute music piece (580.36 seconds, "Reflejos del Sur," Oneness, 2006), both embedded within three resting-state sessions (each lasting 3.3 minutes). The order of the speech and music conditions was counterbalanced across participants. For analysis purposes, the three resting-state sessions were concatenated into a single 10-minute rest condition. During the experiment, sound presentation was facilitated by different setups depending on the location. In the Faraday cage, a Sound Blaster X-Fi Xtreme Audio, a Yamaha P2040 amplifier, and Yamaha loudspeakers (NS 10M) were utilized. Conversely, in the bedroom setting, stimuli were presented using a Sennheiser HD 25 headphone set. Sound stimuli were presented with a 44.1 kHz sample rate and 16-bit resolution, with speech and music excerpts delivered at approximately 75 dBA.

**Preprocessing.** Offline, contact data underwent conversion to virtual channels utilizing a bipolar montage approach, specifically employing a closest-neighbor contact reference. This methodology was chosen to enhance spatial resolution and mitigate passive volume diffusion from neighboring areas (36). Precise localization of the channels was achieved through a procedure akin to the one utilized in the iELVis toolbox (37). Additionally, anatomical localization of the channels was labeled using the Brainnetome Atlas (38). Subsequently, bipolar channels located outside the brain were excluded from the dataset, constituting approximately 3% of the total channels. The remaining data underwent bandpass filtering between 0.1 Hz and 250 Hz. If required, a notch filter was also applied at 50 Hz and harmonics up to 200 Hz to
eliminate power line artifacts. Finally, the data were downsampled to a sampling rate of 500 Hz.

**High-frequency activity (HFa).** The amplitude of high-frequency activity (HFa, 80-120 Hz) was derived by computing the analytic amplitude of four 10-Hz-wide sub-bands spanning the 80 to 120 Hz range, using the Hilbert transform technique. Each sub-band underwent standardization by dividing it by its mean value, and subsequently, all sub-bands were averaged together (39,40). In order to optimize computational efficiency, the data were downsampled to a frequency of 100 Hz.

**Artifact rejection.** Channels exhibiting a variance exceeding 2 times the interquartile range (IQR), a non-parametric estimate of the standard deviation, were identified and labeled as artifacted channels. Subsequently, for the remaining channels, time-segments displaying activity surpassing 5 standard deviations were flagged as artifacted and replaced with the mean of the signal. This approach ensured preservation of the temporal structure of the data, which is essential for evaluating inter-subject correlations.

**Statistics.** To identify significant group differences between speech and music listening and resting state conditions, we employed a cluster-based permutation test (41) to the results obtained using the MDF (Results section) and dFC (Supplementary Material section) metrics. In particular, we performed the test on the differences between paired conditions. For the results obtained using the MDF metric, only vertices with data values exceeding a cluster-forming threshold of 1 were used to form clusters. For the dFC metric, a range of cluster-forming thresholds from 0.5 to 3.0, with a step size of 1, were tested, confirming that no significant differences were found across these varying levels. In each case, we used 10000 permutations. Furthermore, to evaluate potential differences in switching rates among the three conditions (Fig.3D, left), we conducted a Wilcoxon test between the paired distributions corresponding to the three tasks.

**Acknowledgement:** This work was co-funded by the European Union (ERC, SPEEDY, ERC-CoG-101043344).

## References

1. Zatorre R, Belin P, Penhune VB. Structure and function of auditory cortex: music and speech. Trends Cogn Sci (Regul Ed). 2002 Jan 1;6(1):37–46.

2. Krumhansl CL. Cognitive foundations of musical pitch. Oxford University Press; 2001.

3. Carpentier SM, McCulloch AR, Brown TM, Faber SEM, Ritter P, Wang Z, et al. Complexity Matching: Brain Signals Mirror Environment Information Patterns during Music Listening and Reward. J Cogn Neurosci. 2020 Apr;32(4):734–45.

4. te Rietmolen N, Mercier M, Trébuchon A, Morillon B, Schön D. Speech and music recruit frequency-specific distributed and overlapping cortical networks. BioRxiv. 2022 Oct 9;
5. Neri M, Runfola C, te Rietmolen NAG, Sorrentino P, Schon D, Morillon B, et al. Neuronal avalanches in naturalistic speech and music listening. BioRxiv. 2023 Dec 16;

6. Birbaumer N, Lutzenberger W, Rau H, Braun C, Mayer-Kress G. Perception of music and dimensional complexity of brain activity. Int J Bifurcation Chaos. 1996 Feb;06(02):267-78.

7. Bressler SL, Menon V. Large-scale brain networks in cognition: emerging methods and principles. Trends Cogn Sci (Regul Ed). 2010 Jun;14(6):277-90.

8. Hebb DO. The organization of behavior; A neuropsychological theory. Am J Psychol. 1949 Oct;63(4):633.

9. Jirsa V. Structured Flows on Manifolds as guiding concepts in brain science. In: Viol K, Schöller H, Aichhorn W, editors. Selbstorganisation – ein Paradigma für die Humanwissenschaften: Zu Ehren von Günter Schiepek und seiner Forschung zu Komplexität und Dynamik in der Psychologie. Wiesbaden: Springer Fachmedien Wiesbaden; 2020. p. 89-102.

10. Shine JM, Breakspear M, Bell PT, Ehgoetz Martens KA, Shine R, Koyejo O, et al. Human cognition involves the dynamic integration of neural activity and neuromodulatory systems. Nat Neurosci. 2019 Feb;22(2):289-96.

11. Pillai AS, Jirsa VK. Symmetry Breaking in Space-Time Hierarchies Shapes Brain Dynamics and Behavior. Neuron. 2017 Jun 7;94(5):1010-26.

12. Gonzalez-Castillo J, Fernandez IS, Lam KC, Handwerker DA, Pereira F, Bandettini PA. Manifold learning for fMRI time-varying functional connectivity. Front Hum Neurosci. 2023 Jul 11;17:1134012.

13. Nguyen CH, Karavas GK, Artemiadis P. Inferring imagined speech using EEG signals: a new approach using Riemannian manifold features. J Neural Eng. 2018 Feb;15(1):016002.

14. Proix T, Delgado Saa J, Christen A, Martin S, Pasley BN, Knight RT, et al. Imagined speech can be decoded from low- and cross-frequency intracranial EEG features. Nat Commun. 2022 Jan 10;13(1):48.

15. Orepic P, Truccolo W, Halgren E, Cash SS, Giraud A-L, Proix T. Neural manifolds carry reactivation of phonetic representations during semantic processing. BioRxiv. 2024 Jan 21;

16. Saha P, Fels S. Hierarchical Deep Feature Learning for Decoding Imagined Speech from EEG. AAAI. 2019 Jul 17;33(01):10019-20.

17. Ji X, Han J, Jiang X, Hu X, Guo L, Han J, et al. Analysis of music/speech via integration of audio content and functional brain response. Inf Sci (Ny). 2015 Mar;297:271-82.

18. Stephen EP, Li Y, Metzger S, Oganian Y, Chang EF. Latent neural dynamics encode temporal context in speech. Hear Res. 2023 Sep 15;437:108838.

19. Darcourt J. Music's manifold emotions for mental well-being: a topological approach. Res Sq. 2023 Nov 28;

20. Ding N, Patel AD, Chen L, Butler H, Luo C, Poeppel D. Temporal modulations in speech and music. Neurosci Biobehav Rev. 2017 Oct;81(Pt B):181-7.

21. Martin S, Millán JDR, Knight RT, Pasley BN. The use of intracranial recordings to decode human language: Challenges and opportunities. Brain Lang. 2019 Jun;193:73-83.

22. Mercier MR, Dubarry A-S, Tadel F, Avanzini P, Axmacher N, Cellier D, et al. Advances in human intracranial electroencephalography research, guidelines and good practices. Neuroimage. 2022 Oct 15;260:119438.

23. Gardner RJ, Hermansen E, Pachitariu M, Burak Y, Baas NA, Dunn BA, et al. Toroidal topology of population activity in grid cells. Nature. 2022 Feb;602(7895):123-8.

24. Kuo P-C, Chen Y-S, Chen L-F. Manifold decoding for neural representations of face viewpoint
and gaze direction using magnetoencephalographic data. Hum Brain Mapp. 2018 May;39(5):2191-209.

25. Rué-Queralt J, Stevner A, Tagliazucchi E, Laufs H, Kringelbach ML, Deco G, et al. Decoding brain states on the intrinsic manifold of human brain dynamics across wakefulness and sleep. Commun Biol. 2021 Jul 9;4(1):854.

26. Busch EL, Huang J, Benz A, Wallenstein T, Lajoie G, Wolf G, et al. Multi-view manifold learning of human brain-state trajectories. Nat Comput Sci. 2023 Mar 27;3(3):240-53.

27. Chaudhuri D, Li D, Che Y, Shafran E, Gerton JM, Zang L, et al. Enhancing long-range exciton guiding in molecular nanowires by H-aggregation lifetime engineering. Nano Lett. 2011 Feb 9;11(2):488-92.

28. Nieh EH, Schottdorf M, Freeman NW, Low RJ, Lewallen S, Koay SA, et al. Geometry of abstract learned knowledge in the hippocampus. Nature. 2021 Jul 5;595(7865):80-4.

29. Rubin A, Sheintuch L, Brande-Eilat N, Pinchasof O, Rechavi Y, Geva N, et al. Revealing neural correlates of behavior without behavioral measurements. Nat Commun. 2019 Oct 18;10(1):4745.

30. Low RJ, Lewallen S, Aronov D, Nevers R, Tank DW. Probing variability in a cognitive map using manifold inference from neural dynamics. BioRxiv. 2018 Sep 16;

31. Kato S, Kaplan HS, Schrödel T, Skora S, Lindsay TH, Yemini E, et al. Global brain dynamics embed the motor command sequence of Caenorhabditis elegans. Cell. 2015 Oct 22;163(3):656-69.

32. Preti MG, Bolton TA, Van De Ville D. The dynamic functional connectome: State-of-the-art and perspectives. Neuroimage. 2017 Oct 15;160:41-54.

33. Areshenkoff C, Gale DJ, Standage D, Nashed JY, Flanagan JR, Gallivan JP. Neural excursions from manifold structure explain patterns of learning during human sensorimotor adaptation. eLife. 2022 Apr 19;11.

34. Ansmann G. Efficiently and easily integrating differential equations with JiTCODE, JiTCDDE, and JiTCSDE. Chaos. 2018 Apr;28(4):043116.

35. Rackauckas C, Nie Q. Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory. Discrete Continuous Dyn Syst Ser B. 2017;22(7):2731-61.

36. Mercier MR, Bickel S, Megevand P, Groppe DM, Schroeder CE, Mehta AD, et al. Evaluation of cortical local field potential diffusion in stereotactic electro-encephalography recordings: A glimpse on white matter signal. Neuroimage. 2017 Feb 15;147:219-32.

37. Groppe DM, Bickel S, Dykstra AR, Wang X, Mégevand P, Mercier MR, et al. iELVis: An open source MATLAB toolbox for localizing and visualizing human intracranial electrode data. J Neurosci Methods. 2017 Apr 1;281:40-8.

38. Fan L, Li H, Zhuo J, Zhang Y, Wang J, Chen L, et al. The human brainnetome atlas: A new brain atlas based on connectional architecture. Cereb Cortex. 2016 Aug;26(8):3508-26.

39. Ossandón T, Vidal JR, Ciumas C, Jerbi K, Hamamé CM, Dalal SS, et al. Efficient “pop-out” visual search elicits sustained broadband y activity in the dorsal attention network. J Neurosci. 2012 Mar 7;32(10):3414-21.

40. Vidal JR, Freyermuth S, Jerbi K, Hamamé CM, Ossandon T, Bertrand O, et al. Long-distance amplitude correlations in the high y band reveal segregation and integration within the reading network. J Neurosci. 2012 May 9;32(19):6421-34.

41. Maris E, Oostenveld R. Nonparametric statistical testing of EEG- and MEG-data. J Neurosci Methods. 2007 Aug 15;164(1):177-90.
Supplementary Material

<figure><img src="image_6.png" /></figure>

**Results using dynamic Functional Connectivity as a metric.** On the left, the means of the distributions of the standard deviation σ of the dynamic Functional Connectivity (dFC) matrices plotted against different numbers of time window length (in seconds) employed in the computation of the dFC. No significant differences are observed between speech, music, and rest conditions (P>0.07). This supports the effectiveness of our methodology.

<figure><img src="image_7.png" /></figure>

**Results using three independent components to build the MDF.** On the left, the means of the distributions of the standard deviation σ of the MDF matrices plotted against different numbers of time window length (in seconds) employed in the computation of the MDF. Again, speech and music listening display higher complexity than rest, while music listening maintains simpler dynamics compared to speech. This confirms the trend observed when using only two independent components in the construction of the MDF. On the right, the same distributions corresponding to the case in which the window length used to compute the MDF is set to 15 seconds