{'fraction_non_alphanumeric': 0.035050995412639736, 'fraction_numerical': 0.020464971273326506, 'mean_word_length': 4.665783497350492, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Influence of music on empirical brain networks

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 and its increase just before transitions between different parts of the musical form (musical high-level events). The results also suggest a split 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. Finally, we discuss the general modalities of the influence of music on the human brain.


INTRODUCTION
Dealing with the dynamics of neuronal 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). Recently, partial synchronization Schöll (2021) has become a reference point for the explanation of the first-night effect (Tamaki et al., 2016) and unihemispheric sleep (Rattenborg et al., 2000(Rattenborg et al., , 2016Mascetti, 2016;Ramlow et al., 2019). 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 (Joris et al., 1994;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 (ERPs), 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 increase of synchronization represents musical large-scale form perception. In more detail, the authors observed the strongest brain activity at frequencies in the beta-band (about 13 − 30 Hz) and in the gamma-band (about 30 − 120 Hz): In these 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 neuronal dynamics during perception (Bader, 2020). On the other hand, the general influence of sound on empirical brain networks has been investigated recently . 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. 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.
An intriguing synchronization phenomenon in networks is relay (or remote) synchronization between layers which are not directly connected, and interact via an intermediate (relay) layer (Leyva et al., 2018). The simplest realization of such a system is a triplex network where a relay layer in the middle acts as a transmitter between the two outer layers. Remote 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 in the individual layers of a three-layer multiplex network. It has been shown that the three-layer structure of the network allows for (partial) synchronization of chimera states 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 Sect. 2 we discuss the transformation of music to a neuronal input signal using a detailed cochlea model. In Sect. 3 we introduce the neuronal network model based upon empirical connectivities with neuronal input to the auditory cortex generated by music. In Sect. 4 we introduce some methods to characterize the neuronal output. Sect. 5 presents the results of the computer simulations and discusses the dynamical scenarios. Sect. 6 presents a comparison with experiments on human subjects, and Sect. 7 finally concludes.

FROM SOUND TO NEURAL SPIKES
The transformation of sound into neural spikes is the subject of much current research (Bader, 2015(Bader, , 2017(Bader, , 2018. Music, speech, or any sound enters through the outer and middle ear as sound pressure, then acting on the oval window of the cochlear. 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. 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.
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/192000 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 so-called 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 Fig. 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 best-frequency 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 time-delay 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 Fig. 1C.

NEURONAL NETWORK MODEL
We consider an empirical structural brain network as shown in Fig. 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 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  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 Fig. 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 intra-hemispheric coupling compared to inter-hemispheric coupling, is highlighted ( Fig. 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. Figure2B displays the time-series of impulses which are supplied to the brain by means of the auditory cortex. These neuronal impulses were obtained by a method developed by Bader from the hip 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 neuronal 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 seconds and was released in 2002 by American rapper Nas.
hop music song One Mic, composed by the American rapper Nas and released in 2002. Bader's method models the transformation of acoustic mechanical waves to electrical spikes in the cochlea, the part of the inner ear involved in hearing (Bader, 2015(Bader, , 2017(Bader, , 2018. During the transition from 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 neuronal 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: 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 inter-hemispheric 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 activator-inhibitor 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 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 Sect. 2, it has a physical dimension in seconds. On the other hand, the FitzHugh-Nagumo model uses dimensionless time. In order to compare our simulations with real data, 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 (corresponding to dimensionless frequency f FHN ≈ 0.4) 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

METHODS
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 use the length of the input signal I(t) as ∆T . 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, the Kuramoto order parameter is calculated by means 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. 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 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).
The input signal I(t) is obtained from the original music song One Mic by the cochlea model described in Sect. 2 (see Fig. 1). The song has a length of about 4.5 minutes 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 . The local dynamics of our model, the FitzHugh-Nagumo oscillator, has no explicit time scale. Its intrinsic angular frequency is dimensionless and given by ω k = ω FHN = 2πf FHN ≈ 2.51. To include the time signal I(t) correctly in our dimensionless model, we scale the dimensional neuronal input obtained in Sect. 2 by the frequency n b (in Hz), or equivalently, convert the dimensionless output of our simulation into physical dimensions f b according to Eq. (5). This 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).

FREQUENCY BANDS AND COHERENCE
We investigate dynamical scenarios emerging from an external stimulus in the auditory cortices of both hemispheres (k = 41, 86). In order to compare our simulations with the empirical analysis of the influence of music upon the brain Bader, 2014, 2020), 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 Fig. 3. We find a strong non-monotonic 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 Fig. 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). Figure 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 Fig. 4D. It represents a part of the neuronal input signal I(t) constructed from the music song One Mic and shown in Fig. 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 Fig. 3: With increasing value of n b in panels (A)-(C), the time scale of the simulated neuronal 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 Fig. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 4A for n b = 5, might imply a higher value of γ according to Eq. (8), Fig. 4B shows that it is more important that R(t) and I(t) show a similar temporal modulation, as in Fig. 4B for n b = 30. Despite the averaging over 250 simulations over the whole simulation time in Fig. 3, the time segment in Fig. 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.
In the right column of Fig. 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 mean-field 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 Fig. 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 Fig. 4B, the mean phase velocity of the auditory cortex approaches ω having a bigger impact on the dynamics of the whole system than in Fig. 4A for n b = 5 Hz. Remarkable is the fact of a dynamical asymmetry shown by the mean phase velocities in Fig. 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
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;Tallon-Baudry and Bertrand, 1999;Rodriguez et al., 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). 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 − 30 Hz) 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, p.114). 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.
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 Mélen, 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 (EDM), 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, such as 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 that corresponds to the perceived density may also contribute to the emergence of musical form (Bader, 2013;Hartmann and Bader, 2020;Bader, 2021;Lindenbaum et al., 2021;Bader et al., 2021).
Based on the aforementioned correlations between the processes associated with the perception of musical form and neuronal synchronization we expect the dynamics of neural synchronization to correspond to the dynamics of the amplitude of the stimulus, i.e., synchronization to be high when the amplitude of the piece is high and low when the amplitude is low. We would most likely expect this synchronization to occur in this frequency spectrum due to the prominence of the gamma-band, but because synchronization has been observed in other frequency bands as well, we also examine the other frequency bands.
In an experiment, we have recorded the electroencephalogram (EEG) from human scalps to examine the perception of music large-scale form (Fig. 5). 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). 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. For a more detailed description of the experimental procedure, technical details and pre-processing, see (Hartmann and Bader, 2020).
In Fig. 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). The continuous wavelet transform of a 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 may also relate to the fact that all human senses relate physics to perception in a logarithmic way. 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 Fig. 6 and α is the smallest wavelet scale.
For each electrode pair of these nine data sets filtered in this way, the synchronization is calculated in the next step. Thus, we can analyze the synchronization dynamics as a function of the frequency bands. As a synchronization measure the linear cross-correlation, for simplicity taken without time delays, is used, which in fact results in the Pearson correlation coefficient r. 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): where x, y denotes the mean of x, y, respectively. Since we aim to reveal synchronization dynamics on the level of musical form, we calculate the correlation within successive 1-second time windows for each possible pair of electrodes of each wavelet-filtered dataset, which results in 32*31/2*9 = 4464 time series with a resolution of 1 second, and each time series has a length of 270 s corresponding to the stimulus length (see Fig. 7). In recent decades, various methods for measuring synchronization have been introduced (Bastos and Schoffelen, 2015;Blinowska, 2011). 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).
In order to relate this huge number of time series to the amplitude modulation of the stimulus, we first average the amplitude of the stimulus and the correlation coefficients calculated for the 496 electrode pairs and 9 frequency bands within successive 4-second 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 (Fig. 7). In the second step, we correlate all 4464 time series from the synchronization analysis with the amplitude time series of the stimulus. In the third step, we select the 25 time series per frequency band that correlate most strongly with the time series of the amplitude of the stimulus, shown in Fig. 8. Now, we average these 25 time series per frequency band, which results in a single time series of 270 s length for each frequency band, respectively. These averaged time series for each frequency band are correlated over the whole recorded time with the amplitude dynamics of the stimulus (see Fig. 9A). In can be shown that the low and the high gamma-band (FB 2 and FB 3) correlate strongly with the stimulus as expected, but also the slow oscillations (FB 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 over all frequency bands and correlate the resulting time series 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.
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 timeseries as a function of the individual frequency bands is shown in Fig. 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 FB 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 bands FB 6-9 can be interpreted as caused by the interaction of the neocortex with subcortical brain regions. Likewise, the high correlations in bands FB 2-3 are interpreted as activity of the neocortex solely, as expected. The results therefore also suggest a split of musical form-related synchronization between cortical (bands 2-3) and subcortical (bands 6-9) regions.
The high correlations observed in frequency bands 2 and 3 (see Fig. 9A) correspond to a frequency range of 31.25 − 125 Hz (gamma-band). On the other hand in Fig. 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 Fig. 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 Fig. 9B, where the respective frequency bands are averaged from Fig. 3. Both experimental and numerical results show a pronounced maximum of correlation between stimulus and brain dynamics for the gamma-band (bands 2 and 3) in Fig. 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 neuronal 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 A B coherence measure. It has turned out that this coherence measure depends strongly on the frequency band parameter and has its maximum in the corresponding 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 heard. We have experimentally analyzed in detail the influence of real music on the neuronal activity with respect to the common frequency bands in the brain. By means of the Pearson correlation coefficient, 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. Furthermore, the results suggest a split 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. 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 (Kim et al., 2018;Ribeiro et al., 2010;Steyn-Ross and Steyn-Ross, 2010).
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.