Influence of Sound on Empirical Brain Networks

1 Introduction

Synchronization phenomena are well-known regarding dynamical activities of the brain. A high degree of synchronization is related to (slow-wave) sleep [1, 2] or transitions from wakefulness to sleep [3, 4]. Recently, partial synchronization has become a clue to explain the first-night effect [5] and unihemispheric sleep [1, 6–8]. Moreover, synchronized dynamics play an important role in the dynamics of epileptic seizures [9], where the synchronization of a part of the brain causes dangerous consequences for the persons concerned. In contrast, synchronization is also used to explain brain processes which subserve for development of syntax and its perception [10–12]. In general, synchronization theory is highly important to analyze and understand musical acoustics and music psychology [13–17]. While the neurophysiological processes when listening to music remain ongoing research, it is presumed that a certain degree of synchrony can be observed while listening to music and building up expectations. Event-related potentials (ERPs), measured by electroencephalography (EEG) of participants while listening to music, show synchronized dynamics between different brain regions [18, 19]. These studies indicate that the increase of synchronization represents musical large-scale form perception. Moreover, it has been observed that areas of the whole brain are involved regarding neuronal dynamics during perception [10]. Therefore, we propose to investigate the general influence of sound on empirical brain networks. We model the spiking dynamics of the neurons by the paradigmatic FitzHugh–Nagumo model, and investigate possible partial synchronization patterns induced by an external sound source, which is connected to the auditory cortex of the human brain. Furthermore, It is a well-known fact that an important feature of musical sound perception is tonal fusion [20]. Although sound has in general a rich overtone spectrum, subjects perceive only one musical pitch which is a fusion of all partials of the spectrum. Against this background, we concentrate our general study on an external sound source with an amplitude and a single frequency, neglecting the complexity of music and its distinct effects in different frequency bands within the brain oscillations. Within the scope of this work, we have restricted ourselves to a minimal model with no node-specific behavior to reveal the impact of a periodic perturbation.

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 [21]. 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 [22–26]. 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 [27–31]. 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 will be affected by the distance between regions and therefore a stimulation 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 [32]. 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 [33]. 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 [18, 19]. Therefore, we defer the consideration of time delays for now.

2 Model

We consider an empirical structural brain network shown in Figure 1 where every region of interest is modeled by a single FitzHugh–Nagumo (FHN) oscillator.

FIGURE 1. (color online) 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 [36], but re-labeled such that $k=1,\dots ,45$ and $k=46,\dots ,90$ correspond to the left and right hemisphere, respectively. After [9].

The weighted adjacency matrix $A=\left\{A_{kj}\right\}$ of size 90 × 90, with node indices $k\in N=\left\{\mathrm{1,2},\dots ,90\right\}$ 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 [34], for previous utilization of the structural networks to analyze chimera states see [7, 9, 35]. 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 Diffusion Tensor Imaging (DTI) and subsequently divided into 90 predefined regions according to the Automated Anatomical Labeling (AAL) atlas [36]. 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 1 are rearranged such that $k\in N_L=\left\{\mathrm{1,2},\dots ,45\right\}$ corresponds to left and $k\in N_R=\left\{46,\dots ,90\right\}$ to the right hemisphere. Thereby the hemispheric structure of the brain, i.e., stronger intra-hemispheric coupling compared to inter-hemispheric coupling, is highlighted (Figure 1).

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 [37]. 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 [38].

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=\mathrm{41,86}$ (temporal sup L/R). The auditory cortex takes part in the spectrotemporal analysis of the inputs passed on from the ear.

Each node corresponding to a brain region is modeled by the FitzHugh–Nagumo (FHN) model with external stimulus, a paradigmatic model for neuronal spiking [39–41]. 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 [42, 43]. Thus the dynamics of the network reads:

where $\epsilon =0.05$ describes the timescale separation between the fast activator variable (neuron membrane potential) $u$ and the slow inhibitor (recovery variable) $v$ [40]. Depending on the threshold parameter $a$, the FHN model may exhibit excitable behavior $\left(\left|a\right|>1\right)$ or self-sustained oscillations $\left(\left|a\right|<1\right)$. 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 external stimulus is modeled by a trigonometric function with frequency ω and amplitude γ and is applied to the brain areas $k=\mathrm{41,86}$ associated with the auditory cortex, i.e. $C_k=1$ if $k=41$ or 86 and zero otherwise. The coupling between the single regions is given by the coupling strength $\sigma$. As we are looking for partial synchronization patterns we fix $\sigma =0.6$ similar to numerical studies of synchronization phenomena during unihemispheric sleep [7] and epileptic seizures [9] where partial synchronization patterns have been observed. The interaction scheme between nodes is characterized by a rotational coupling matrix:

with coupling phase $\varphi =\frac{\pi }{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 [44] as it reduces the stability of the completely synchronized state. Also in the modeling of epileptic-seizure-related synchronization phenomena [9], 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 [45], and gives rise to partial synchronization patterns which are not found otherwise.

3 Methods

We explore the dynamical behavior by calculating the mean phase velocity ${\omega }_k=2\pi M_k/\Delta T$ for each node $k$, where $\Delta T$ denotes the time interval during which M complete rotations are realized. Throughout the paper we use $\Delta T=\mathrm{10,000}$. For all simulations we use initial conditions randomly distributed on the circle $u_k^2+v_k^2=4$. In case of an uncoupled system (σ = 0), the mean phase velocity (or natural frequency) of each node is ${\omega }_k={\omega }_{FHN}\approx 2.6$. Furthermore we introduce hemispheric measures that characterize the degree of synchronization of the sub-networks and give complementary information. First, the spatially averaged mean phase velocity is:

Thus $\overline{\omega }$ 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 ${\theta }_k$ that can be obtained from the standard geometric phase ${\tilde{\varphi }}_k\left(t\right)=\arctan \left(v_k/u_k\right)$ by a transformation which yields constant phase velocity ${\dot{\theta }}_k$. For an uncoupled FHN oscillator the function $t\left({\tilde{\varphi }}_k\right)$ is calculated numerically, assigning a value of time $0<t\left({\tilde{\varphi }}_k\right)<T$ for every value of the geometric phase, where T is the oscillation period. The dynamical phase is then defined as ${\theta }_k=2\pi t\left({\tilde{\varphi }}_k\right)/T$, which yields ${\dot{\theta }}_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. Additionally, we calculate the temporal mean of the Kuramoto order parameter

to estimate the general dynamical behavior of the system over time. Similarly, the temporal mean $\langle \Omega \left(t\right)\rangle$ of the collective frequency Ω of the mean field [46], defined by

can be considered, and compared with the spatially averaged mean phase velocity.

4 Synchronization Regions

We investigate synchronization scenarios emerging from an external periodic stimulus in the auditory cortices of both hemispheres ($k=\mathrm{41,86}$). Figure 2 shows synchronization scenarios of an empirical structural brain network in dependence of the frequency $\omega$ and amplitude γ of the external stimulus. The light colored regions in Figure 2A indicates synchronized dynamics, whereas the darker colors indicate desynchronized dynamics. There is a light colored stripe for ω = 2.6 which indicates a Kuramoto order parameter $\langle R\rangle \approx 0.8$ and a light colored tongue starting at $\omega =2.4,\gamma =0.04$. The hatched region in Figure 2A stands for a low standard deviation < 0.1 of the temporal mean of the Kuramoto order parameter $\langle R\rangle$. It indicates the absence of strong fluctuations of $R\left(t\right)$ and therefore a constant high level of synchrony in time. Figure 2B shows the drop of the spatially averaged mean phase velocity $\overline{\omega }$ in case of coherent dynamics in the synchronization regions of Figure 2A. In the upper region, $\overline{\omega }$ takes over the value of the frequency $\omega$ of the external stimulus, whereas in the synchronization tongue $\overline{\omega }$ keeps its value of $\overline{\omega }=2.4$.

FIGURE 2. (color online) Synchronization tongues in brain network with external stimulus: (A) The temporal mean of the Kuramoto order parameter $\langle R\rangle$ for simulation time $\text{Δ}T=\mathrm{10,000}$ and (B) the spatially averaged mean phase velocity $\overline{\omega }$ in the parameter plane of the frequency ω of the external stimulus and its amplitude γ. The light color in panel (A) stands for synchronization and the darker color for desynchronization. In the hatched region the standard deviation of $\langle R\rangle$ is less than 0.1, which indicates the absence of strong fluctuations of R in time. The dynamics of the four marked dots in each panel are shown in Figures 3A–D, 4A–D. Other parameters are given by σ = 0.6, $\mathit{ϵ}=0.05$, a = 0.5, and $\varphi =\frac{\pi }{2}-0.1$.

It turns out that by taking the frequency $\omega$ of the external stimulus as a control parameter, one can change the level of synchrony of the system. Figure 3 depicts the details of the transition to synchronization for increasing values of the frequency ω of the external stimulus. Fixing the amplitude $\gamma =0.06$, we take a closer look on the temporal evolution of R and the mean phase velocities in the system for different regions in Figure 2: In Figure 3A the temporal evolution of the Kuramoto order parameter is similar to the system behavior without external stimulus, i.e., it exhibits large temporal fluctuations. In the right column the phase velocities of all nodes are plotted, the horizontal grey dotted line indicates the temporal average of the collective mean-field frequency Ω. Only the phase velocity of the auditory cortex follows the frequency of the external driving stimulus ω = 2.3 and therefore is lower than the frequency of the other nodes ${\omega }_k\approx 2.8$. Increasing the external frequency to ω = 2.4 yields an abrupt transition to a synchronized state. In Figure 3B the Kuramoto order parameter $R\approx 0.95$ and the mean phase velocities indicate a synchronous dynamical behavior, which agrees with the collective frequency Ω of the mean-field (grey dotted horizontal line). With a further increment to ω = 2.5, the system loses synchrony (see Figure 3C) and enters the region between the two synchronization regions in Figure 2A. For ω = 2.6 in Figure 3D, which corresponds to the natural frequency of the uncoupled oscillators, the system regains synchronization, though the Kuramoto order parameter with $R\approx 0.8$ is lower than in the synchronization tongue. Remarkable is the fact of a dynamical asymmetry shown by the mean phase velocities. While the nodes of the right hemisphere exhibit an equal mean phase velocity, a part of the left hemisphere exhibits a faster dynamic similar to dynamics of unihemispheric sleep studied in [7]. In such states one hemisphere is synchronized, whereas the other hemisphere is partly desynchronized.

FIGURE 3. (color online) Dynamical scenarios: dynamics inside and outside the synchronization regions (marked as black dots in Figure 2) by the Kuramoto order parameter R (left column) and the mean phase velocities ${\omega }_k$ (right column) for increasing values of the frequency $\omega$ of the external stimulus ω = 2.30 (A), ω = 2.44 (B), ω = 2.50 (C), and ω = 2.60 (D) for fixed amplitude $\gamma =0.06$. The vertical dashed line in the right column separates the left and right hemisphere; the horizontal grey dotted line indicates the temporal average of the mean-field frequency ω. The red dots mark the nodes of the auditory cortical regions (k = 41, 86). Other parameters are as in Figure 2.

For a better insight, Figure 4 shows the space-time plot of the variable u_k for the corresponding parameter values in Figure 3. In Figures 4B,D, the dynamics inside the two synchronization regions is depicted. The perturbation in the mean phase velocity profile in the right panel of Figure 3D, can be detected also in the corresponding perturbations in Figure 4D. Comparing Figures 4A,C, we can see an increase of synchronized time segments. This increase will be analyzed quantitatively in more detail in the inset of Figure 5.

FIGURE 4. (color online) Synchronized and desynchronized dynamics: Shown are the space-time plot of the variable $u_k$ inside and outside the synchronization regions (marked as black dots in Figure 2) for increasing values of the frequency ω of the external stimulus ω = 2.30 (A), ω = 2.44 (B), ω = 2.50 (C), and ω = 2.60 (D) for fixed amplitude $\gamma =0.06$. The panels correspond to the panels in Figure 3. Other parameters are as in Figure 2.

FIGURE 5. (color online) Transition scenarios: (A) temporal mean of the Kuramoto order parameter $\langle R\rangle$ (dark blue) and the spatially averaged mean phase velocities $\overline{\omega }$ (light orange) in dependence on the frequency $\omega$ of the external stimulus for a fixed amplitude $\gamma =0.052$. The vertical bars indicate the standard deviation of the temporal mean of the Kuramoto order parameter and the spatially averaged mean phase velocities, respectively. As input nodes, the auditory cortices k = 41, 86 are chosen. In case of a different input (k = 1, 45) the corresponding light grey curves are shown in panel (A). The inset in panel (A) depicts ${\rho }_s=\frac{N_s}{\text{Δ}{T}_L}$, the number $N_s$ of synchronized time intervals $\left(R\left(t\right)>0.8\forall t\right)$ divided by a simulation time of $\text{Δ}{T}_L=\mathrm{30,000}$ for values of the frequency ω between the two synchronization regions. The vertical bars denote the standard deviation of the length of these synchronized time intervals. (B)$\langle R\rangle$ for a larger range of driving frequencies $\omega$, showing higher resonance tonges. Other parameters are as in Figure 2.

5 Transition to Synchronization

There are two frequencies which play an important role for the dynamics of the system. On the one hand, in Figure 2A a broad synchronization region is located at a frequency $\omega \approx 2.6$, which is the frequency of the uncoupled FHN oscillator ${\omega }_{\text{FHN}}$. Although the external stimulus effects only the two auditory nodes (k = 41,86), we can observe a transition to synchronization of the whole system approaching $\omega \approx 2.6$ already for small values of the amplitude $\gamma >0.004$. On the other hand, we can detect a synchronization tongue with a lower boundary at $\omega \approx 2.4$ and an upper boundary increasing linearly with the amplitude γ. In contrast to the first, smooth transition, we can find here a sharp transition to synchronized dynamics, similar to a first order transition, depicted by the high contrast of the boarder of the synchronization tongue in Figure 2A. In this synchronization tongue, the nodes oscillate with an equal mean phase velocity (see Figure 3B), but there are phase differences between them, as indicated by $0.95<R\left(t\right)<1$ and shown in the phase-time plot in Figure 4B. Using the fact that $u_j$/$v_j$ and $u_k$/$v_k$ are on the same limit cycle in the phase space and have the same mean phase velocity, the phase differences in the coupling term of Eq. 1 can be effectively summed up in following way:

where $\Delta t_{eff }\ll 1$ denotes the effective sum of the time intervals of all phase differences. Neglecting $\cos \varphi \ll 1$ and setting $\sin \varphi \approx 1$, Eq. 1 reads for $k\ne \mathrm{41,86}$:

The local dynamics of Eq. 1 is governed by a slow-fast system (FitzHugh-Nagumo oscillator), where the slow part essentially determines the period of the oscillations. Hence, considering the slow motion on the falling branches of the u-nullcline (${\dot{u}}_k=0$) by inserting the second equation into the first one

the time derivative of the falling branches yields with ${\dot{v}}_k$ from Eq. 8

The separation of the variables gives

where $dt$ can be integrated over one oscillation period T. As shown in [47], this leads in case of synchronization to a linear dependence of the oscillation period $T_{sync}={\displaystyle {\int }_0^T\text{d}t}$ on the effective sum of the phase differences proportional to $\Delta t_{eff }$. For incoherent distribution of the phases of each node $k$ (see Figure 4D), the phase differences between the single nodes are also strongly distributed and thus $\Delta t_{eff }\approx 0$. In this case, the natural frequency of the uncoupled system plays an important role, provided that the mean phase velocity of all oscillators is still almost equal as in case of Figure 2D.

This could explain on one side the fact that we observe a synchronization tongue at $\omega \approx 2.4$ (which is smaller than the frequency of an uncoupled oscillator ${\omega }_{FHN}\approx 2.6$), and on the other side, the linear boundaries of the synchronization tongue for increasing amplitude γ. The increase of γ yields an increase of the sum of the phase differences in the coupling term of Eq.1 and therefore an increase of the effective sum of the time intervals $\Delta t_{eff }$.

In Figure 5A, both transitions are depicted in dependence on the frequency ω for a fixed amplitude $\gamma =0.052$. We can see an abrupt increase and decrease of the temporal mean of the Kuramoto order parameter $\langle R\rangle$ before and after $\omega \approx 2.4$, respectively. In contrast, in approaching the upper synchronization region starting from $\omega \approx 2.6$, $\langle R\rangle$ increases more slowly than at the transition to the synchronization tongue ($\omega \approx 2.4$). In case of synchronization the standard deviation of $\langle R\rangle$, displayed by the vertical bars, is smaller than in case of desynchronized dynamics. That holds also for the spatially averaged mean phase velocities $\overline{\omega }$, which in case of synchronization takes over the lower value of the frequency ω of the external stimulus. Also for $\omega >2.6$, $\overline{\omega }$ is equal to $\omega$, whereas the standard deviation of $\overline{\omega }$ increases linearly with $\omega$. In contrast, there is no effect on the system for $\omega <2.4$. Neither $\langle R\rangle$ nor $\overline{\omega }$ show a different behavior for such values of $\omega$. The high value of the standard deviation of $\langle R\rangle$ stands for dynamics as shown in Figure 3A, where the Kuramoto order parameter $R\left(t\right)$ is fluctuating over its whole bandwidth $R\in \left[\mathrm{0,1}\right]$. Simulations show that for ω > 3.0 the dynamical behavior of the system becomes similar to that with ω = 2.3. For both parameter intervals of ω, there is no effect on the system. Simulations show also that a similar transition to synchronization at ω = 2.6 can be found for higher harmonics, i.e., multiple values of ω = 2.6. In Figure 5B, we can identify synchronization regions for $\omega =5.2,7.8$, and 10.4 becoming less pronounced for increasing $\omega$, i.e., having a smaller extension in the plane of ω and γ. In contrast, we could not detect repeated synchronization tongues of ω for multiple values of ω = 2.4. This indicates the existence of two different synchronization mechanisms.

The existence of two synchronization regions depends on the choice to which nodes the external stimulus is supplied. In case of a different input, for instance k = 1,45 in contrast to k = 41,86, the light grey curves in Figure 5A depict the corresponding dependence of the Kuramoto order parameter $\langle R\rangle$ and the spatially averaged mean phase velocities $\overline{\omega }$ upon the frequency ω of the external stimulus. The synchronization region at $\omega \approx 2.4$ is missing here and only one synchronization region remains ($\omega >2.5$).

In the following, we analyze the region between the two synchronization areas in more detail. Figure 6 depicts the dynamical behavior when we approach the synchronization region by increasing the frequency ω of the external stimulus in the neighborhood of the synchronization region at ω = 2.6. For ω = 2.47 in Figure 6A, the time series of the Kuramoto order parameter shows familiar temporal fluctuations with only short episodes of synchronization $\left(R\left(t\right)>0.8\right)$. In [9] the authors define the threshold $R=0.8$ as the onset of an epileptic seizure. By increasing the frequency ω, one can increase the quantity of these episodes, as well as their duration. Figure 6D with ω = 2.51 features much longer duration of synchronized episodes, moreover the duration of the single episodes are comparable in length. This transition in Figures 6A–D can be also seen in Figure 5A. The inset of Figure 5A confirms the increasing regularity between the two synchronization regions by depicting ${\rho }_s=\frac{N_s}{\Delta T_L}$ vs. $\omega$, where $N_s$ is the number of synchronized time intervals ($R\left(t\right)>0.8\forall t$) and $\Delta T_L=\mathrm{30,000}$ is the simulation time. The vertical bars denote the standard deviation of the length of these synchronized time intervals. With increasing $\omega$ not only the number of synchronized time intervals is increasing, but the standard deviation of their duration is decreasing. For $\omega$ = 2.6 we enter the synchronization region, where the value of ${\rho }_s$ drops due to the nearly consistently synchronized dynamics.

FIGURE 6. (color online) Transition scenarios: The increase of the regularity and duration of synchronized time intervals shown by the temporal evolution of the Kuramoto order parameter R (left column) and the mean phase velocities ${\omega }_k$ (right column) for increasing values of the frequency ω of the external stimulus ω = 2.47 (A), ω = 2.49 (B), ω = 2.51 (C), and ω = 2.53 (D) for fixed amplitude $\gamma =0.06$. The vertical dashed line in the right column separates the left and right hemisphere, and the horizontal grey dotted line indicates the temporal average of the mean-field frequency Ω. The red colored dots indicate the nodes of the auditory cortical regions (k = 41,86). Other parameters as in Figure 2.

Finally, the mean phase velocities in the right column of Figure 6 display the transition to frequency synchronization. While the frequency of the two driven nodes (k = 41, 86) converges to the frequency of an uncoupled FHN oscillator ${\omega }_{FHN}\approx 2.6$, also the frequencies of all the other nodes are adjusted, especially those with a much higher frequency (k = 18, 63).

6 Conclusion

We have investigated the influence of an external sound source on the dynamics of a network with empirical structural connectivity. It has been found that depending on the frequency and amplitude of the sound source, synchronization can be induced in the dynamics of the system. We have shown that two frequencies play an important role for synchronization, particularly the natural frequency of the uncoupled oscillator and the frequency of the coupled system. Moreover, the degree of synchronization is gradually increased when the frequency of the uncoupled oscillator or multiple values of it are approached. Furthermore, we have analyzed the linear dependence of the synchronization borders upon the amplitude of the external sound, which can also be characterized as the volume of the sound. This has resulted in the observation that the synchronization region can be enlarged by increasing volume. We have demonstrated the dynamical behavior of the system in the transition to synchronization. By tuning the frequency of the external sound appropriately, we have shown that the level of synchrony can be increased.

These results are in accordance with experiments of Bader’s group [18, 19] that music induces a certain degree of synchrony in the human brain. This group has shown that listening to music can have remarkable influence on the brain dynamics, in particular, a periodic alternation between synchronization and desynchronization. Moreover, such an alternation 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 [45], exhibiting hysteresis and avalanche phenomena as seen in critical phenomena and phase transitions [48–50]. By choosing appropriate parameters, we have reported an intriguing dynamical behavior regarding the transition to synchronization, and have observed the induced alternation between high and low degrees of synchronization. To sum up, an external sound source connected to the brain allows for synchronization dynamics, which may be used to model the effect of music on the human brain.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Author Contributions

JS did the numerical simulations and the theoretical analysis. ES supervised the study. All authors designed the study and contributed to the preparation of the article. All the authors have read and approved the final article.

Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project No. 429685422).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Antonín Škoch and Jaroslav Hlinka (National Institute of Mental Health, Klecany, Czech Republic) for providing the sample structural connectivity matrices, and to Rolf Bader and Lenz Hartmann (University of Hamburg) for stimulating discussions.

References