Modelling the perception of music in brain network dynamics

We analyze the influence of music in a network of FitzHugh-Nagumo oscillators with empirical structural connectivity measured in healthy human subjects. We report an increase of coherence between the global dynamics in our network and the input signal induced by a specific music song. We show that the level of coherence depends crucially on the frequency band. We compare our results with experimental data, which also describe global neural synchronization between different brain regions in the gamma-band range in a time-dependent manner correlated with musical large-scale form, showing increased synchronization just before transitions between different parts in a musical piece (musical high-level events). The results also suggest a separation in musical form-related brain synchronization between high brain frequencies, associated with neocortical activity, and low frequencies in the range of dance movements, associated with interactivity between cortical and subcortical regions.


Introduction
Dealing with the dynamics of neural networks, one repeatedly encounters the phenomenon of synchronization. In the brain, a high degree of synchronization is related to (slow-wave) sleep (Steriade et al., 1993;Rattenborg et al., 2000) or transitions from wakefulness to sleep (Schwartz and Roth, 2008;Moroni et al., 2012). Often, only a part of the brain is synchronized. This phenomenon of so-called partial synchronization Schöll (2021) has recently become a reference point for the explanation of unihemispheric sleep (Rattenborg et al., 2000(Rattenborg et al., , 2016Mascetti, 2016;Ramlow et al., 2019) and the first-night effect (Tamaki et al., 2016), which describes troubled sleep in a novel environment. Furthermore, synchronized dynamics plays an integral role in the dynamics of epileptic seizures (Gerster et al., 2020), where the synchronization of a part of the brain causes dangerous consequences for the persons concerned. By contrast, synchronization is also used to explain brain processes serving the development of syntax and its perception (Koelsch et al., 2013;Large et al., 2015;Bader, 2020). Generally, synchronization theory is of great importance for the analysis and understanding of musical acoustics and music psychology (Bader, 2013;Sawicki et al., 2018a;Hou et al., 2020;Shainline, 2020).
Although the neurophysiological processes involved in listening to music are still being researched, it is believed that some degree of synchrony can be observed in listening to music and building expectations. Event-related potentials, measured by electroencephalography (EEG) of participants while listening to music, show synchronized dynamics between different brain regions Bader, 2014, 2020). These studies indicate that the synchronization dynamics represents musical large-scale form perception. The coupling of oscillatory neural signals within the usual frequency bands has been thought to be a mechanism that is related to a broad range of perceptual, sensorimotor, and cognitive processes, such as Gestalt perception and binding (Gray and Singer, 1989;Tallon et al., 1995;Keil et al., 1999;Rodriguez et al., 1999;Tallon-Baudry and Bertrand, 1999;, timing and expectation Meck, 2005, 2009), attention (Womelsdorf and Fries, 2007;Fries, 2009;Nikolić et al., 2013), consciousness (Baars, 2006;Dehaene et al., 2011;Engel and Fries, 2016;Owen and Guta, 2019), or motor functions (Thaut et al., 2015) as well as in music perception (Bhattacharya et al., 2001;Zanto et al., 2005;Bonetti et al., 2021).
According to (Engel and Fries, 2016), oscillatory brain activity is usually clustered into several frequency bands: delta (0.5-3.5 Hz), theta (4-7 Hz), alpha (8-12 Hz), beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) and gamma ( > 30 Hz). Since the gamma-band is the 'youngest' frequency band which has become of interest (from about the late 1990s), the ranges and definitions vary from source to source. Here, we refer to the classification of (Freeman and Quian Quiroga, 2013), who speak of a low gamma range for frequencies above 30 Hz up to 60 Hz, and high gamma for frequencies above 60 Hz up to about 120 Hz. For everything above 120 Hz, we use the term 'fast oscillations' as employed by Buzsáki (2006). The gamma-band frequency range is of particular interest in the context of large-scale synchronization since it is thought to be a mechanism that integrates information from different parts of the cortex. In more detail, for specific frequency bands the increase and decrease of synchronization are following the large-scale form of the listened music in a coherent way. Moreover, it has been observed that areas of the whole brain are involved in neural dynamics during perception (Bader, 2020).
The musical form as the hierarchically highest level of musical structure and its perception is related to some of the mentioned processes above (Lerdahl and Jackendoff, 1990;Hartmann and Bader, 2020). Perceptually, notes, bars, and phrases are grouped and integrated into a high-level part of the form by the Gestalt laws (Leman, 1997;Deutsch, 2013;Neuhaus, 2013;Deliége and Melen, 2014). The contrast of the form's parts, such as the concatenation of verse and chorus in a song, the sonata form of classical music, or the continuous night-long tension build-up and decay in Techno, House or Electronic Dance Music, characterize the musical form and the learned knowledge about the underlying structures leads to the build-up of expectation and their fulfillment as well as to modulated attention. On an emotional level, this can be expressed in the terms of tension and relaxation (Koelsch, 2014;Lehne and Koelsch, 2015). Also, the transition from "potential energy" (expectations) into "kinetic energy" (dancing) as proposed by (Kurth, 1931) can be related to the processing of musical form in the sense of entrainment of neurons in the motor cortex by neurons from the auditory cortex (Thaut et al., 2015).
The characteristic of contrasting parts can be revealed not only by music analysis using pen and paper but also by different computational methods by the music information retrieval discipline, like the amplitude of a piece of music that corresponds to the subjective perception of loudness. Also other properties of the stimulus, such as the spectral centroid that corresponds to the perceived brightness of a sound, or the fractal correlation dimension (Grassberger and Procaccia, 1983a,b) corresponding to the perceived density and thereby representing the complexity of a piece of music, are drivers of the musical form (Bader, 2013;Hartmann and Bader, 2020;Bader, 2021;Bader et al., 2021;Linke et al., 2021).
Recently, the general influence of sound on a dynamical system with complex network connectivities (derived from empirical Diffusion Tensor Imaging (DTI) measurements) has been investigated . It has been shown that an external sound source, which is connected to the auditory cortex of the human brain, induces partial synchronization patterns. Nevertheless, this study has neglected the complexity of music and its distinct effects in different frequency bands within the brain oscillations. There are a variety of recognized modeling approaches with respect to neural systems in general (Kacprzyk and Pedrycz, 2015;Bassett and Sporns, 2017;Bassett et al., 2018;Petkoski et al., 2018;Petkoski and Jirsa, 2019) and related to music in particular (Friston and Friston, 2013). In this paper, we model the spiking dynamics of the neurons by the paradigmatic FitzHugh-Nagumo model, and investigate possible coherence between the dynamics of the brain network and an external music source, which is connected to the auditory cortex of the human brain. Moreover, we present experimental data which we successfully reproduce numerically with the help of our network model, which combines simple node dynamics with complex network connectivities derived from empirical measurements.
An intriguing synchronization phenomenon in multilayer networks is relay synchronization between layers which are not directly connected, and interact via an intermediate (relay) layer (Leyva et al., 2018). Multilayer networks can give a general framework to describe and model real life examples of various systems, e.g., the two hemispheres of the brain or two cortical regions connected by the hippocampus (Gollo et al., 2011). Relay synchronization, a regime where pairs of nodes synchronize despite their large distances on the network graph, has been shown to depend on the network symmetries (Bergner et al., 2012;Nicosia et al., 2013;Gambuzza et al., 2013;Zhang et al., 2017a,b). Recently the notion of relay synchronization has been extended from completely synchronized states to partial synchronization patterns. It has been shown that the multilayer structure of a network allows for (partial) synchronization in the outer layers via the relay layer (Sawicki et al., 2018b,c;Sawicki, 2019;Winkler et al., 2019;Drauschke et al., 2020;. Going towards more realistic models, time-delay plays an important role in the modeling of the dynamics of complex networks. In brain networks, the communication speed is affected by the distance between regions and therefore a stimulus applied to one region needs time to reach a different region. In such delayed system, it is possible to predict if the effects of stimulation remain focal or spread globally (Muldoon et al., 2016). More generally, time delays due to propagation over the white-matter tracts have been shown to organize the brain network synchronization dynamics for different types of oscillatory nodes (Petkoski and Jirsa, 2019). Within the scope of this paper, we focus on the requirements for a simple model to exhibit partial synchronization patterns, which have been experimentally observed Bader, 2014, 2020). Therefore, we defer the consideration of time delays for now.
This article is organized as follows. In Section 2, we discuss the transformation of music to a neural input signal using a detailed cochlea model. In Section 3, we introduce the neural network model based upon empirical connectivities with neural input to the auditory cortex generated by music. In Section 4, we introduce some methods to characterize the neural output. Section 5 presents the results of the computer simulations and discusses the dynamical scenarios. Section 6 presents a comparison with experiments on human subjects, and Section 7 finally concludes.

From sound to neural spikes
The transformation of sound into neural spikes is the subject of much current research (Tritsch et al., 2010;Mizrahi et al., 2014;Bader, 2015Bader, , 2017Bader, , 2018Guo et al., 2021). Music, speech, or any sound enters through the outer and middle ear as sound pressure, then acting on the oval window of the cochlea. The movement of the oval window is then transferred to a pressure in the lymph liquid of the cochlea surrounding the basilar membrane, which again acts on the basilar membrane, causing traveling waves there. Due to spatial differences in stiffness and damping on the membrane, sinusoidal waves with a single frequency show an increase in amplitude up to a point with maximum amplitude, the position of the so-called best-frequency, with a fast decay afterwards. Therefore, different positions on the basilar membrane represent different frequencies, making the cochlea a Fourier analyzer. The stereocilia on the basilar membrane at the position of respective best-frequency are then transferring the mechanical energy into neural spikes. The frequency distribution on the basilar membrane is logarithmic. Movements of neighboring frequencies lead to interactions, causing roughness perception up to a frequency band of a musical major third. These bands are called critical bands, and the basilar membrane consists of 24 such bands. The spikes leaving the respective bands are fed into the auditory pathway, consisting of several neural nuclei, where the nucleus cochlearis or the trapezoid body are the first two. The interaction between these neural nuclei is manifold with several feedback loops and binaural connections (Schofield, 2011) ending at the auditory cortex on both hemispheres. Still up to the A1 region of the auditory cortex, the critical bands are maintained, where neural connections of higher nuclei are connected to bands on the basilar membrane, which is called tonotopy.
Many auditory features are present, extracted, or perceived already in this pathway, like sound localization, pitch, or timbre (Lyon and Shamma, 1996), although research has not concluded on further processing in the cortex (Bader, 2021). Music perception of larger temporal content, like song or sonata form, are not part of processing in the auditory pathway up to the cortex, as far as we know. Still the feedback loops within the pathway are both directions, up and down, afferent and efferent, so e.g. there is one connection down from the cortex to the cochlea with only one nucleus in between, tuning the basilar membrane tension through efferent nerves, according to cortex activity (Schofield, 2011).
Up to now, no model of the whole auditory pathway exists on a detailed neural level. The model used in this paper therefore concentrates on main findings, i.e., the transition from sound to neural spikes, the tonotopy of neural connections up to the cortex, as well as partial synchronization of phases in the cochlea, which are also present as coincidence detection in the auditory pathway. A Finite-Difference Time Domain (FDTD) physical model of the cochlea is used (Bader, 2015). The basilar membrane is about 3.5 cm long and only between 0.1-0.12 cm wide, so it is more a rod than a membrane. Therefore, the present model assumes a differential equation of a membrane like with basilar membrane displacement u along a one-dimensional axis x, basilar membrane stiffness K(x) = 2 × 10 9 e −3.4x dyn/cm 3 changing along x, and linear mass density μ(x) = m/A(x) with mass m over cross section A again changing along the basilar membrane and A(x) = 0.1 cm × (0.1 cm + 0.02 cm × x/l) with basilar membrane length l = 3.5 cm taking into account the slight widening of the basilar membrane over its length. The boundary conditions of the basilar membrane are homogeneous Dirichlet boundary conditions which do not allow for displacements on Frontiers in Network Physiology frontiersin.org the boundaries, but any derivative is allowed in accordance with the physiological conditions. Comparison between a membrane and a rod model shows no considerable differences, therefore a rod model is used. Here d is damping, and f(t) is the driving force of the lymph fluid which drives the basilar membrane.
To calculate the spikes omitted by the cochlea, the recording of the musical piece used is fed into the cochlea model. Here the amplitudes of the digital musical sound file are taken as sound pressures acting on the oval window of the cochlea and therefore immediately on the peri-and endolymph around the basilar membrane. As the speed of sound in the lymph (~1,500 m/s) is much larger than the speed of waves on the basilar membrane which is between~100 m/s at the oval window and down to1 0 m/s at the helicotrema, an instantaneous action of the pressure at the oval window on the basilar membrane is reasonable and known as long-wave approximation (de Boer, 1991). This holds for frequencies up to~4 kHz, where pitch perception stops and humans only hear a very high sound. This approximation is used in the model. It leads to the force f(t) in Eq. 1 which represents the amplitudes of the digital musical sound file acting instantaneously on all points of the basilar membrane at each time point respectively. It is interesting to see that the traveling wave on the basilar membrane is therefore not caused by an external input slowly traveling through the cochlea but by the intrinsic solution of the inhomogeneous differential equation of the basilar membrane driven by a periodic force over its whole length instantaneously.
Depending on the brain region, neurological measurements reveal different time scales (Spitmaan et al., 2020). In our work we choose 50 ms as a time integration step as this is consistent with a characteristic time scale in music as well as in visual perception. In music 50 ms correspond to the second integration time, below which two events cannot be distinguished one from another. This leads to a threshold of 20 Hz, above which musical pitches are perceived and below which adjacent events are heard as rhythms. In vision, 18-24 frames per second lead to a continuous visual perception, again corresponding to about 50 ms time intervals. Therefore, in terms of hearing and seeing, the brain seems to update perceptional input on a time-scale of 50 ms (Bader, 2013).
The transition between mechanical displacement and electrical spike is performed using two conditions according to literature (Hubbard and Mountain, 1996). A neural spike at one point X on the basilar membrane at time τ is excited if two conditions hold.
Condition (2a) means a maximum shearing of two nervous fibers as a necessary condition to an opening of the ion channels at the fibers. This only happens with a positive slope, as only then the stereocilia are driven away from each other. With a negative slope the cilia are getting closer and therefore no stress appears at the tip links between them. This corresponds to the rectification process in gammatone filter banks. Condition (2b) is a temporal maximum positive peak of the basilar membrane displacement. It is the temporal equivalent to the spatial condition of a maximum acceleration, where the tip link between the cell and its neighboring cells is most active.
To calculate the spikes omitted by the cochlea, the recording of the musical piece used is fed into the cochlea model. Therefore, the original piece, available as a digital recording of 44.1 kHz sample rate (CD-Quality) is upsampled to 192 kHz to meet Finite-Difference Time Domain (FDTD) stability criteria. The cochlea model is then run with a time-discretization step of Δt = 1/192,000 s. Each time when a neural spike appears, the time point, strength, and critical band of the spike is stored. Therefore, after processing, a time series I(t) of all spikes leaving the cochlea is obtained. Figure 1A displays an example of an artificially generated socalled tone complex with f 0 = 475 Hz and ten partial tones (harmonics) with amplitudes 1/m where m = 1, 2, 3, . . . , 10. The respective spike output of the basilar membrane model is shown in Figure 1B. Each time when the sound wave has a maximum amplitude, a pressure pulse is traveling over the basilar membrane, which emits electrical spikes at respective bestfrequency positions on the membrane in accordance with the frequencies in the activating sound. As traveling waves on the membrane start at the basal end, next to the oval window, where high frequencies have their best-frequency location, and travel down the membrane towards the upper end, the helicotrema, where low frequencies are located, low frequencies show a timedelay with respect to higher frequencies. If the spikes of all critical bands are summed up for a certain point in time, a time series I(t) of all neural spikes leaving the cochlea can be generated, as exemplarily shown in Figure 1C. The simplification that the output of the cochlea model is summed up at one time point is motivated by the results of (Joris et al., 1994): In an experiment with cats, the authors could show that the scattered output of the cochlea is synchronized in the trapezoid body.

Neural network model
In this section, we introduce an empirical structural brain network as shown in Figure 2A where every region of interest is modeled by a single FitzHugh-Nagumo (FHN) oscillator. The weighted adjacency matrix A = {A kj } of size 90 × 90, with node indices k, j ∈ N = {1, 2, . . . , 90} was obtained from averaged diffusion-weighted magnetic resonance imaging data measured in 20 healthy human subjects. For details of the measurement procedure including acquisition parameters, see (Melicher et al., 2015), for previous utilization of the structural networks to analyze chimera states see (Chouzouris et al., 2018;Ramlow et al., 2019;Gerster et al., 2020;Schöll, 2021). The data were analyzed using probabilistic tractography as implemented in the Frontiers in Network Physiology frontiersin.org 04 FMRIB Software Library, where FMRIB stands for Functional Magnetic Resonance Imaging of the Brain (www.fmrib.ox.ac.uk/ fsl/). The anatomic network of the cortex and subcortex is measured using Diffusion Tensor Imaging (DTI) and subsequently divided into 90 predefined regions according to the Automated Anatomical Labeling (AAL) Atlas (Tzourio-Mazoyer et al., 2002), see Table 1. Each node of the network corresponds to a brain region. Note that in contrast to the original AAL indexing, where sequential indices correspond to homologous brain regions, the indices in Figure 2A are rearranged such that k ∈ N L = {1, 2, . . . , 45} corresponds to left and k ∈ N R = {46, . . . , 90} to the right hemisphere. Thereby the hemispheric structure of the brain, i.e., stronger intrahemispheric coupling compared to inter-hemispheric coupling, is highlighted ( Figure 2A).
The structural connectivity matrices serve as a realistic input for modeling, rather than as exact information concerning the existence and strength of each connection in the human brain. The pipeline for constructing such connectivity information using diffusion tractography is known to face a range of challenges (Schilling et al., 2019). While some estimates of the strength and direction of structural connections from measurements of brain activity can in principle be attempted, the relation of these can vary dramatically with (experimentally unknown) parameters of the local dynamics and coupling function (Hlinka and Coombes, 2012).
The auditory cortex is the part of the temporal lobe that processes auditory information in humans. It is a part of the auditory system, performing basic and higher functions in hearing and is located bilaterally, roughly at the upper sides of the temporal lobes, i.e., corresponding to the AAL indexing k = 41, 86 (temporal sup L/R). The auditory cortex takes part in the spectrotemporal analysis of the input passed on from the ear. Figure 2B displays the time-series of impulses which are supplied to the brain by means of the auditory cortex. These neural impulses were obtained by the method of Bader described in Section 2 (Bader, 2015(Bader, , 2017(Bader, , 2018. Here, in contrast to Figure 1, a real piece of music was used, namely the hip hop music song One Mic, composed by the American rapper Nas and released in 2002. During the transition from acoustic mechanical to electrical excitation within the cochlea, synchronization appears to improve perception of pitch, speech, or localization. The sampling rate of these impulses obtained by Bader's method is f s = 192 kHz. Each node corresponding to a brain region is modeled by the FitzHugh-Nagumo (FHN) model with external stimulus, a paradigmatic model for neural spiking (FitzHugh, 1961;Nagumo et al., 1962;Bassett et al., 2018). Note that while the FitzHugh-Nagumo model is a simplified model of a single neuron, it is also often used as a generic model for excitable media on a coarse-grained level (Chernihovskyi et al., 2005;Chernihovskyi and Lehnertz, 2007). Thus the dynamics of the network reads:  (Bader, 2015), where the vertical axis represents the cochlea position with best-frequency f in Hz indicated, i.e., categorized into 24 so-called critical bands. (C) Time series I(t) of the sum of all spike weights leaving the cochlea at a certain time t. Note that the first 5 ms are transients.
With k ∈ N H where N H denotes either the set of nodes k belonging to the left (N L ) or the right (N R ) hemisphere. Parameter ϵ = 0.05 describes the timescale separation between the fast activator variable (neuron membrane potential) u and the slow inhibitor (recovery variable) v (FitzHugh, 1961). Depending on the threshold parameter a, the FHN model may exhibit excitable behavior (|a| > 1) or self-sustained oscillations (|a| < 1). We use the FHN model in the oscillatory regime and thus fix the threshold parameter at a = 0.5 sufficiently far from the Hopf bifurcation point. The coupling within the hemispheres is given by the coupling strength σ while the coupling between the hemispheres is given by the interhemispheric coupling strength ς. As we are looking for partial synchronization patterns we fix σ = 0.7 and ς = 0.15 similar to numerical studies of synchronization phenomena during unihemispheric sleep (Ramlow et al., 2019) where partial synchronization patterns have been observed. The interaction scheme between nodes is characterized by a rotational coupling matrix: with coupling phase ϕ π 2 − 0.1, causing primarily an activatorinhibitor cross-coupling. This particular scheme was shown to be crucial for the occurrence of partial synchronization patterns in ring topologies (Omelchenko et al., 2013) as it reduces the stability of the completely synchronized state. Also in the FIGURE 2 (A) Model for the hemispheric brain structure: Weighted adjacency matrix A kj of the averaged empirical structural brain network derived from twenty healthy human subjects by averaging over the coupling between two brain regions k and j. The brain regions k, j are taken from the Automated Anatomic Labeling Atlas (Tzourio-Mazoyer et al., 2002), but re-labeled such that k = 1, . . . , 45 and k = 46, . . . , 90 correspond to the left and right hemisphere, respectively. After (Gerster et al., 2020). (B) Time-series of the neural input signal I(t) obtained from the music song One Mic transformed by a method developed by Bader (Bader, 2020). The song has a length of about 270 s and was released in 2002 by American rapper Nas.

Frontiers in Network Physiology
frontiersin.org modeling of epileptic-seizure-related synchronization phenomena (Gerster et al., 2020), where a part of the brain synchronizes, it turned out that such a cross-coupling is important. The subtle interplay of excitatory and inhibitory interaction is typical of the critical state at the edge of different dynamical regimes in which the brain operates (Massobrio et al., 2015;Shi et al., 2022), and gives rise to partial synchronization patterns which are not found otherwise.
The external stimulus I(t) describes the impulses evoked by the music piece One Mic by Nas and is applied to the brain areas k = 41, 86 associated with the auditory cortex, i.e., C k = 1 if k = 41 or 86 and zero otherwise. Since I(t) is a time series which is calculated from a real piece of music, see Section 2, it has a physical dimension in seconds. On the other hand, the FitzHugh-Nagumo model has no explicit time scale. Its intrinsic angular frequency is dimensionless and given by ω k = ω FHN = 2πf FHN ≈ 2.51 (corresponding to dimensionless frequency f FHN ≈ 0.4). In order to compare our simulations with real data and include the time signal I(t) correctly in our dimensionless model, we must transform the dimensionless time units of the FHN oscillator model to real time units by comparing the FHN oscillation period of a single FHN oscillator T ≈ 2.5 to the characteristic frequencies n b in Hz of an empirical time series. Depending upon the frequency band n b (in Hz) chosen, the simulation time is converted to real time by 1 s = 2.5n b simulation time units, or the simulated frequency (in Hz) is In this way, the parameter n b effectively removes the time scale from the input, but on the other hand, it can also be seen as creating a link between our dimensionless model and the input signal I(t).

Synchrony measures
We explore the dynamical behavior by calculating the mean phase velocity ω k = 2πM k /ΔT for each node k, where ΔT denotes the time interval during which M k complete rotations are realized. Throughout the paper, we denote the length of the input signal I(t) as ΔT. For the numerical integration an adaptive Runge-Kutta integration method has been applied (python scipy: solve_ivp, RK45). For all simulations we use initial conditions randomly distributed on the circle u 2 k + v 2 k 4 and a transient time of t trans = 10,000 before the input signal I(t) is supplied to the system. In case of an uncoupled system (σ = ς = 0), the mean phase velocity (or natural frequency) of each node is ω k = ω FHN = 2πf FHN ≈ 2.51.
First, we introduce the spatially averaged mean phase velocity: Thus ω corresponds to the mean phase velocity averaged over the left and right hemisphere.
Second, we take advantage of an abstract dynamical phase θ k that can be obtained from the standard geometric phaseφ k (t) arctan(v k /u k ) by a transformation which yields constant phase velocity _ θ k . For an uncoupled FHN oscillator the function t(φ k ) is calculated numerically, assigning a value of time 0 < t(φ k ) < T for every value of the geometric phase, where T is the oscillation period. The dynamical phase is then defined as θ k 2πt(φ k )/T, which yields _ θ k const. Thereby identical, uncoupled oscillators have a constant phase relation with respect to the dynamical phase. By means of the dynamical phase θ k we can calculate the Kuramoto order parameter where the fluctuations of the order parameter R caused by the FHN model's slow-fast time scales are suppressed and a change in R indeed reflects a change in the degree of synchronization. The Kuramoto order parameter may vary between 0 and 1, where R = 1 corresponds to complete phase synchronization, and small values characterize spatially desynchronized states. Third, we introduce a new measure which specifies the coherence between the Kuramoto order parameter and the input signal by using the time average of the Kuramoto order parameter weighted with the input signal to quantify the overlap of coherent episodes (R large) with large input signals, averaged over time. The coherence γ is maximum if the synchronization is large whenever the signal is large. It is small if the overall synchronization is low, or if the modulation of the synchronization in time is not in phase with the modulation of the input signal amplitude. For γ = 0 the Kuramoto order parameter and the input signal do not overlap at any time point. An increased value of γ ∈ [0, 1] means increased overlap between the Kuramoto order parameter and the input signal. The motivation for introducing the measure γ lies in the fact that in the human brain the increase and decrease of synchronization follows the large-scale form of the listened music in a coherent way Bader, 2014, 2020). Fourth, we make use of the Pearson correlation coefficient r, a linear cross-correlation, for simplicity taken without time delays. This is widely used as a non-directed measure of the strength of the correlation between two variables or sequences {x 1 , x 2 , . . . , x n } and {y 1 , y 2 , . . . , y n } (Glantz, 2002;Bastos and Schoffelen, 2015;Guevara Erra et al., 2017): where x, y denotes the mean of x, y, respectively. In recent decades, various methods for measuring synchronization have been introduced (Blinowska, 2011;Bastos and Schoffelen, 2015). The advantage of the Pearson correlation coefficient r is that it allows for easy and efficient calculation of the linear correlations between two variables or time series, and the results are very similar to those obtained by other common methods such as the phase-locking value . For a comparison of the different synchronization measures see (Jalili et al., 2014).
The input signal I(t) is obtained from the original music song One Mic by the cochlea model described in Section 2 (see Figure 1). The song has a length of about 4.5 min and the sampling rate of the obtained input signal is given by f s = 192 kHz. Sampling is the reduction of a continuous-time signal to a discrete-time signal, e.g., the conversion of a sound wave (a continuous signal) to a sequence of samples (a discrete-time signal). The sampling rate f s is then the average number of samples obtained in one second. According to the Nyquist criterion, the frequency information of I(t) is then band-limited to f b < 1 2 f s .

Frequency bands and coherence
Next, we investigate dynamical scenarios emerging from an external stimulus in the auditory cortices of both hemispheres (k = Frontiers in Network Physiology frontiersin.org 41, 86). In order to compare our simulations with the empirical analysis of the influence of music upon the brain Bader, 2014, Hartmann andBader, 2020, see also Section 6), we may choose different frequency bands n b , and hence a different scaling of the time in the external stimulus. This can be visualized by plotting the coherence measure γ in dependence on the characteristic frequency n b (in Hz), see Figure 3. We find a strong nonmonotonic behavior of γ(n b ) and it turns out that by taking the frequency band n b of the external stimulus as a control parameter, one can change the level of coherence between the system dynamics and the external stimulus. Although the standard deviation of the coherence measure is relatively large for an ensemble size of 200 simulations (indicated by the vertical bars), we find a pronounced maximum of the coherence γ for n b = 12-48 Hz corresponding to the gamma-band of brain waves (f b ≈ 30-120 Hz) shown in Figure 3 by purple shading. This means that for that frequency n b the level of synchronization follows the external signal most closely. It is in agreement with what has been observed in empirical brain analysis of the perception of music Bader, 2014, 2020).
Figures 4A-C depicts the details of the change of the time series of the Kuramoto order parameter R(t) with increasing values of the frequency band n b of the external stimulus I(t), which is shown in Figure 4D. It represents a part of the neural input signal I(t) constructed from the music song One Mic and shown in Figure 2B. We take a closer look at the temporal evolution of R and the mean phase velocities ω k in the system for different values of n b chosen from three different regimes in Figure 3: With increasing value of n b in panels (A)-(C), the time scale of the simulated neural output in Hz changes from lower to higher frequencies f b which is also seen in the temporal fluctuations of R(t). Furthermore we observe on the one hand an increasing amplitude of the temporal fluctuations of R. On the other hand, the temporal average of the Kuramoto order parameter R decreases with increasing n b , marked by a horizontal grey dotted line in the left column: While for a small value of n b = 5 Hz in Figure 4A the Kuramoto order parameter R assumes rather large values, and small values R < 0.2 are not reached, for high values of n b = 90 Hz in Figure 4C rather small values of R are measured. This trend can be seen by means of the temporal average of the Kuramoto order parameter R. For n b = 30 Hz in Figure 4B, the temporal average of R takes a value ≈ 0.5 and the time evolution shows regular oscillations between low (R < 0.2) and high values (R > 0.8). This aspect will be further discussed in the next section, since it can also be observed in experiments.
As shown in Figure 3, in the case of n b = 30 Hz the coherence γ is maximum. Even though a higher value of the temporal average of R(t), as observed in Figure 4A for n b = 5, might imply a higher value of γ according to Eq. (8), Figure 4B shows that it is more important that R(t) and I(t) show a similar temporal modulation, as in Figure  4B for n b = 30. Despite the averaging over 250 simulations over the whole simulation time in Figure 3, the time segment in Figure 4B shows such a similarity in the modulation: We can see simultaneous drops of R(t) < 0.1 and I(t) < 0.1 for example at t ≈ 138, 140, 150, whereas the values in between are higher, even if they fluctuate.  Figure 2B. The vertical dashed line in the right column separates the left and right brain hemisphere; the red dots mark the nodes of the auditory regions (k = 41, 86). The horizontal grey dotted line indicates the temporal average of the Kuramoto order parameter R in the left column, and the spatial average of the mean-field frequency ω in the right column. Other parameters are as in Figure 3.

Frontiers in Network Physiology frontiersin.org
In the right column of Figure 4 the dimensionless mean phase velocities ω k of all nodes are plotted, the horizontal grey dotted line indicates the spatial average, i.e., the collective meanfield frequency ω, which does not change for different n b since it is determined by the intrinsic collective dynamics. In contrast, the node dynamics of the auditory regions (k = 41, 86), indicated by red dots, depends on n b since it receives the external input signal which has a higher frequency in dimensionless units if the time is scaled in larger units 1/n b . For n b = 5 Hz in Figure 4A, the mean phase velocity of the auditory cortex is higher compared to the spatial average of the collective mean-field frequency ω. For n b = 30 Hz in Figure 4B, the mean phase velocity of the auditory cortex approaches ω having a bigger impact on the dynamics of the whole system than in Figure 4A for n b = 5 Hz.
Remarkable is the fact of a dynamical asymmetry shown by the mean phase velocities in Figure 4C: While the nodes of the right hemisphere exhibit equal mean phase velocity, i.e., they are frequency synchronized, the left hemisphere remains desynchronized and exhibits on average faster dynamics. This may indicate that regardless of the input I(t) the system can exhibit partial synchronization. Such behavior is similar to the dynamics of unihemispheric sleep studied in (Ramlow et al., 2019), where no external input has been applied to the dynamical system. In such states one hemisphere is synchronized, whereas the other hemisphere is partially desynchronized.

Comparison with experiments
Based on the correlations between the processes associated with the perception of musical form and neural synchronization, we expect the dynamics of neural synchronization to correspond to the amplitude dynamics of the stimulus. Again, the musical amplitude corresponds to perceived loudness, and is calculated as integration of energy over time intervals. Then synchronization between different brain regions is high when the amplitude of the musical piece is high, and synchronization is low when the amplitude of the piece is low. We expect such brain synchronization to be strong due to the prominence of the gamma-band in perception of musical parameters.
In an experiment, we have recorded the electroencephalogram (EEG) from human scalps to examine the perception of music large-scale form (see Figure 5) 1 . 25 musically skilled subjects listened to the song One Mic from the artist Nas three times each. The song was released in 2001 on his Album Stillmatic on Columbia Records. The electroencephalogram (EEG) signals were recorded with a sample rate of 500 Hz from 32 electrodes, positioned following the 10-20 method of placement (Jasper, 1958). In this experiment, we are focused on the temporal dynamics of synchronization related to the time span of the musical form and therefore do not take advantage of methods for the inverse modeling of EEG data (Schoffelen and Gross, 2009;Palva and Palva, 2012).
After artifact correction, recorded data for each channel has been averaged over subjects and trials to obtain a grand average of 75 trials for each channel to increase the signal to noise ratio and enhance event-related potentials. This type

FIGURE 5
Recorded and averaged electroencephalogram (EEG) data: top and middle plot show recorded EEG time series after pre-processing for one electrode (Fp1) from two different participants. The bottom plot shows the time series of the same electrode averaged over 25 subjects and three trials.
1 We have taken into account the usual guidelines regarding ethical procedure (informed consent). The subjects were mainly found and recruited through the Institute of Systematic Musicology Hamburg and had instrumental lessons on at least one instrument (mean duration 10.0 years, standard deviation 4.6 years) or corresponding experience as DJ. They participated in accordance with local ethics committee guidelines.
Frontiers in Network Physiology frontiersin.org of averaging reveals evoked potentials (in contrast to induced potentials) and is related to the presented stimulus in a classical event-related potential manner (Tallon-Baudry et al., 1996;Tallon-Baudry and Bertrand, 1999;Zanto et al., 2005). We are aware that our choice for evoked potentials pushes subjective, individual brain activity that is not stimuluslocked into the background. Indeed, it was found that this subjective, individual brain activity, often referred to as 'noise', contains valuable information that is lost when averaging over many subjects (Tallon-Baudry and Bertrand, 1999). On the other side, recent studies on this issue have shown strong overlap between subjects' brain activity (Hasson et al., 2004;Dmochowski et al., 2012;Abrams et al., 2013;Kaneshiro et al., 2021). Therefore, we choose to take advantage of the improvement in the signal-to-noise ratio over the disadvantage of the individual portion of the perception. Individual perception might be subject to future studies. Also note that the choice of using a correlation analysis between single electrodes is not including redundant synchrony due to overlap of electrical fields between electrodes, since the positions of the electrodes do not differ over measurement time. Therefore, the differences in correlation strength between different electrodes cannot be explained by spurious synchrony (Holsheimer and Feenstra, 1977;Kayser and Tenke, 2006;Bhavsar et al., 2018). For a more detailed description of the experimental procedure, technical details and pre-processing, see (Hartmann and Bader, 2020).
In Figure 6, all channels have been decomposed into nine independent frequency bands that correspond approximately to the frequency bands mentioned above by using a continuous wavelet transformation with a Mexican Hat wavelet (Freeman and Quian Quiroga, 2013). In contrast to a bandpass filter with a subsequent Hilbert transform, using a Mexican hat wavelet for filtering is fast and efficient since one can decompose the recorded EEG data into the desired frequency bands in one step by defining the number of octaves. The continuous wavelet transform of a uniformly sampled sequence {x 1 , x 2 , . . . x n } = {x (t 0 ), x (t 0 + Δt), . . ., x (t 0 + (n − 1)Δt)} is given by where s ∈ R corresponds to the frequency of the EEG band and u = 1, . . . , n labels the wavelet coefficients with the number n of analyzed sample points defining the time window of observation. As wavelet function ψ a Mexican Hat wavelet is used, given by where σ is the width of the wavelet. The EEG bands used align very well with a musical scale, where each higher band doubles the frequency of its respective lower band, corresponding to a musical octave. Please note that this relation might only be at chance, still it Frontiers in Network Physiology frontiersin.org may also relate to the fact that all human senses relate physics to perception in a logarithmic way (Schneider, 2018). It is therefore convenient to scale s in the wavelet transform in the same mathematical way as an equal-tempered musical scale like s oct = α 2 oct−1 , where oct ∈ {1, 2, . . . , 9} is the octave number related to the nine frequency bands shown in Figure 6 and α is the smallest wavelet scale. For each electrode pair of these nine data sets filtered in this way, the synchronization is calculated by means of the Pearson correlation coefficient r (see Eq. 9) in the next step. Thus, we can analyze the synchronization dynamics as a function of the frequency bands. Since we aim to reveal synchronization dynamics on the level of musical form, we calculate the correlation within successive 1-s time windows for each possible pair of electrodes of each wavelet-filtered dataset, which results in 32*31/2*9 = 4,464 time series of correlation coefficients representing the synchronization dynamics between electrode-pairs with a resolution of 1 s, and each of these time series has a length of 270 s corresponding to the stimulus length (see Figure 7).
In order to relate this huge number of time series of correlation coefficients to the amplitude dynamics of the stimulus, we first average the amplitude of the stimulus and the correlation coefficients calculated for the 496 electrode pairs and nine frequency bands within successive 4-s windows to avoid minor amplitude fluctuations and obtain a scaling corresponding to about two musical bars that fits to changes related to the musical form (Figure 7). In the second step, we correlate all 4,464 time series of correlation coefficients with the amplitude dynamics of the stimulus. In the third step, we select the 25 time series of correlation coefficients per frequency band that correlate most strongly with the amplitude dynamics of the stimulus, shown in Figure 8. Now, we average these 25 time series of correlation coefficients per frequency band, which results in a single time series of 270 s length for each frequency band, respectively. These averaged time series of correlation coefficients, representing the synchronization dynamics for each frequency band, are correlated over the whole recorded time with the amplitude dynamics of the Example of the synchronization dynamics between two electrodes. Dashed black line: Time series of the Pearson correlation coefficient r calculated for successive 1-s time windows (n = 500 in Eq. 9 between averaged EEG recordings of electrode Fp1 (lower plot in Figure 5) and electrode T7. Blue line: Pearson correlation coefficient averaged over four consecutive 1-s time windows of the dashed black line.

FIGURE 8
Comparison of whole brain synchronization dynamics and representation of the musical form of the stimulus. The black line shows the amplitude dynamics of the stimulus as a representation of the musical form, averaged over each of four consecutive seconds. The blue line shows the average of the 25 correlation time series between two electrodes from each frequency band that correlates most strongly with the amplitude dynamics of the stimulus.

Frontiers in Network Physiology
frontiersin.org stimulus (see Figure 9A). It can be shown that the low and the high gamma-band (frequency bands 2-3) correlate strongly with the stimulus as expected, but also the slow oscillations (frequency bands 7-9) correlate very well (see discussion below). By this, we can reveal how good the synchronization dynamics in each frequency band corresponds to the amplitude dynamics of the stimulus on the level of musical form. In the next step, we average these time series representing the synchronization dynamics for each frequency band and correlate the resulting time series, representing the synchronization dynamics of the whole brain, with the amplitude dynamics of the stimulus as well. These two time series correlate with a Pearson coefficient of r = 0.76. Therefore, we can conclude that the higher the amplitude of the stimulus, the higher the synchronization between the most correlated time series of the different frequency bands. According to (Cohen, 1992), this is a strong effect. As shown in Figure 8, the increased synchrony is not constant during music listening, but rather synchronization dynamics follows the sound amplitude. Note that the correlation between sound amplitude (perceived loudness) or other parameters like brightness or fractal correlation dimension (see inset of Figure 9A) and brain synchronization is not trivial. First, brain synchronization appears at frequencies much lower than most musical frequencies. Secondly, synchronization appears with multiple perceptual parameters. Thirdly, increasing, e.g., the sound amplitude might lead to an increase of the network amplitude, but here it leads to an enhanced synchronization, pointing to a highly nonlinear process in the network, caused by the activity of the brain when perceiving sound.
It is interesting to note that the correlation with the stimulus is highest when the time series from all frequency bands are averaged. The correlation coefficient of the averages of the 25 most correlated time-series as a function of the individual frequency bands is shown in Figure 9A. It shows two regimes of high correlation, separated by a frequency band (FB 5) with low correlation. Here, the central nervous system in the spinal cord and its relation to the locomotor system are expected to be responsible for the dynamics in the frequency bands 6-9 due to their frequency range close to walking and dancing (van Noorden and Moelants, 1999). Note that the electroencephalogram (EEG) recordings are performed on the skull, and therefore represent the brain dynamics of the neocortex which is interacting with the brain stem. Therefore, the high correlations between synchronization and musical form in frequency bands 6-9 can be interpreted as caused by the interaction of the neocortex with subcortical brain regions. Likewise, the high correlations in frequency bands 2-3 are interpreted as activity of the neocortex solely, as expected. The results therefore also suggest a separation of musical form-related synchronization between cortical (frequency bands 2-3) and subcortical (frequency bands 6-9) regions.
The high correlations observed in frequency bands 2-3 for the sound amplitude (see Figure 9A) as well as for the fractal correlation dimension (see inset of Figure 9A) correspond to a frequency range of 31.25-125 Hz (gamma-band). On the other hand in Figure 3, the strongest coherence between the Kuramoto order parameter (measure for global neural synchronization) and the external input can be found for n b = 10-40 Hz. Taking into account that the natural frequency of each node is f FHN ≈ 0.4, we can calculate the corresponding frequency band f b = n b /f FHN . As shown by the upper x-axis in Figure 3, the strongest coherence in our model can be observed for a frequency band of f b = 40-100 Hz, which agrees with the gamma-band in the brain. For comparison with the experiment, we show the corresponding numerically simulated results in Figure 9B, where the respective frequency bands are averaged from Figure 3. Both experimental and numerical results show a  (Grassberger and Procaccia, 1983a,b) has been used for the calculation of r. (B) Numerically simulated coherence γ between network dynamics and external stimulus, where the corresponding frequency bands are averaged from Figure 3. As in Figure 3, the purple shaded regions in both panels indicate the gamma-band (f b ≈ 30-120 Hz), respectively.

Frontiers in Network Physiology
frontiersin.org pronounced maximum of correlation between stimulus and brain dynamics for the gamma-band (frequency bands 2-3) in Figure 9. Note that the second maximum in the experimental data (panel A), which is due to the interaction of the neocortex with subcortical brain regions as discussed above, is absent in the simulated data (panel B) since the computer simulation is only performed for the neocortex, using a cochlea input, but neglecting brain stem activity.

Conclusion
We have investigated the influence of music in a simulated network of FitzHugh-Nagumo oscillators with empirical structural connectivity obtained from healthy human subjects, and have compared it to measured electroencephalogram (EEG) data. We report an increase of coherence between the global dynamics and the input signal induced by a specific music song. We have shown that the level of coherence depends on the frequency band. We have compared our results with experimental data, which describe global neural synchronization between different brain regions in the gamma-band range and its increase just before transitions between different parts of the musical form (musical high-level events). Such synchronization increases before musical large-scale form boundaries, and decreases afterwards, therefore represents musical large-scale form perception.
The transformation of sound into neural spikes takes place in the cochlea, a part of the human ear which is directly connected to the auditory cortex. By means of the basilar membrane, the brain is able to perceive different frequencies organized in so-called critical bands. We have applied a cochlea model to transform a specific music song into an input signal representing neural spikes evoked by the music song. This input signal has then been supplied to a simulated network of neural oscillators with empirical structural connectivity. By the transformation of the dimensionless time units of the oscillator model to real time units, we have investigated dynamical scenarios in dependence on the introduced frequency band parameter. To quantify moreover the overlap between input signal and network dynamics, we have introduced a coherence measure. It has turned out that this coherence measure depends sensitively on the frequency band and has its maximum in the gamma-band. Therefore, depending on the frequency band, coherence can be induced between the dynamics of the system and its input signal.
These results are in accordance with our own and previous experiments Bader, 2014, 2020) where music has also been found to induce a certain degree of synchrony in the human brain. We have shown that listening to music can have a remarkable influence on the brain dynamics, in particular, a periodic alternation between synchronization and desynchronization which is strongly related to the music perceived. We have experimentally analyzed in detail the influence of real music on the neural activity with respect to the common frequency bands in the brain. By means of the Pearson correlation coefficient of the sound amplitude as well as the fractal correlation dimension, we have found the gamma-band to be important for musical form perception. Just as in the computer simulation, we have found a pronounced maximum for this frequency range. Moreover as in simulation, the increased gamma-band synchrony is not constant during music listening in our experiment, but rather synchronization dynamics follows the musical largescale form represented by a perceptual related characteristic of the stimulus, i.e., the amplitude and fractal correlation dimension. Even though we chose a specific piece of music in this study, we expect future work to show that these results can be generalized.
Furthermore, the results suggest a separation in musical formrelated brain synchronization between high brain frequencies, associated with neocortical activity, and low frequencies in the range of dance movements, associated with interactivity between cortical and subcortical regions. Besides, an alternation between synchronization and desynchronization reflects the variability of the system; this can be seen as a critical state between a fully synchronized and a desynchronized state. It is known that the brain is operating in a critical state at the edge of different dynamical regimes (Massobrio et al., 2015;Shi et al., 2022), exhibiting hysteresis and avalanche phenomena as seen in critical phenomena and phase transitions (Ribeiro et al., 2010;Steyn-Ross and Steyn-Ross, 2010;Kim et al., 2018).
By choosing appropriate parameters and measures, we have reported an intriguing dynamical behavior in dependence on the frequency bands, and have observed the induced increase of coherence both in numerical and experimental setups. To sum up, music supplied to the brain allows for a high coherence and correlation between musical input and brain dynamics especially in the gamma-band. This insight may be used to fathom the general modalities of the influence of music on the human brain.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Ethics statement
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. The patients/participants provided their written informed consent to participate in this study.

Author contributions
JS did the numerical simulations and the theoretical analysis, LH has performed the experiments. RB and ES Frontiers in Network Physiology frontiersin.org supervised the study. All authors designed the study and contributed to the preparation of the manuscript. All the authors have read and approved the final manuscript.

Funding
This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project No. 429685422) and the Open Access Publication Fund of TU Berlin.