Partial synchronization in empirical brain networks as a model for unihemispheric sleep

A well-known phenomenon in nature is unihemispheric slow-wave sleep, exhibited by aquatic mammals including whales, dolphins and seals, and multiple bird species. Unihemispheric sleep, as the name suggests, is the remarkable ability to engage in deep (slow-wave) sleep with a single hemisphere of the brain while the other hemisphere remains awake [1–3]. Interestingly, sleep and wakefulness are characterized by a high and low degree of synchronization, respectively [1]. Sleep is associated with specific synchronized oscillations, i.e., sleep spindles and slow oscillations in the thalamocortical system [4]. In addition, arousal- and sleep-promoting neural assemblies undergo collective activity resulting in secretion of sleep-regulating neurotransmitters [5]. While the synchronization processes can differ between adults and children [6], transitions from wakefulness to sleep are widely accompanied by synchronization phenomena [7].

In the human brain the first-night effect, which describes troubled sleep in a novel environment, has been recently related to asymmetric dynamics, i.e., a manifestation of one hemisphere of the brain being more vigilant than the other [8]. Sleep is a dynamical macrostate of the brain that is observed over a wide range of animal species. Sleep is accompanied by a loss of consciousness and conscious perceptions, and muscle activity is reduced or absent. Sleep alternates between rapid-eye-movement (REM) and non-REM stages N1, N2, N3, where the latter are dominated by slow oscillations (1 Hz and below) which can also emerge locally [9,10]. Sleep stage switching dynamics includes wake/sleep asymmetric stochasticity [11], but obeys an underlying control by regulatory circuits forming bistable biological flip-flop switches [12–15], and sleep regulation is coupled to the sleep oscillations of the thalamocortical system [16]. While most animals follow a similar qualitative sleep pattern and fall into sleep with both hemispheres, in certain bird and mammal species sleep can be unihemispheric [3]. It has been speculated that unihemispheric sleep is related to the spontaneous symmetry-breaking phenomenon of chimera states in oscillator networks [17,18]; those states combine spatially coexisting domains of synchronized and desynchronized dynamics [19–23].

While the neurophysiological processes that ensure the existence of this dynamical state of unihemispheric sleep remain largely unknown, it is presumed that a certain degree of structural interhemispheric separation is a necessary condition for this pattern to persist. Therefore, we propose to model unihemispheric sleep by a two-community network of the two hemispheres where the inter-hemispheric coupling strength is smaller than the intra-hemispheric coupling. We model the spiking dynamics of the neurons by the paradigmatic FitzHugh-Nagumo model, and investigate possible partial synchronization patterns.

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

Zoom In Zoom Out Reset image size

The brain network was obtained from diffusion-weighted magnetic resonance imaging data measured in healthy human subjects as part of a larger study focusing on connectivity changes in schizophrenia. For details of the measurement procedure including acquisition parameters, see [24], for previous utilization of the structural networks to analyze chimera states see [25]. 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 [26]. Each node of the network corresponds to a brain region. Indirect information of the white-matter fibers connecting different brain regions is provided by diffusion-weighted Magnetic Resonance Imaging (dMRI) measuring the preferred diffusion direction in each voxel of the brain. Probabilistic tractography then provides for each voxel a set of $n_s = 5000$ streamlines, simulating the possible white-matter fiber tracts. A coefficient P_kj giving the connectivity probability from the k-th to the j-th region is introduced by the proportion of streamlines connecting voxels in region k to voxels of region j on the condition that they originate in region k. Thus, a weighted adjacency matrix of size $90 \times 90$, with node indices $k \in N = \{1,2,\ldots,90\}$ is constructed. Finally, each entry in this adjacency matrix, i.e., the connectivity between every two regions, is averaged over 20 subjects (mean age 33 years, standard deviation 5.7 years, 10 females, 2 left-handed) yielding the average empirical structural brain network $\mathbf{A} = \{A_{kj}\}$. The pipeline for constructing the structural network has been adopted from a previous study of differences in connectivity patterns between healthy subjects and schizophrenia patients [27]. Obtaining such connectivity information using diffusion tractography is known to face a range of challenges [28,29]. Note that in contrast to the original AAL indexing, where sequential indices correspond to homologous brain regions, the indices in fig. 1 are rearranged such that $k \in N_L = \{1, 2, \ldots,45\}$ corresponds to the left hemisphere and $k \in N_R = \{46, \ldots,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. 1(a)). Note that there is a very slight structural asymmetry of the two brain hemispheres. This structural asymmetry is well described in the neuroscience literature and is related to the known asymmetries in localization of psychological functions, such as the prevalence of language functions in the left brain hemisphere in humans. Moreover, this lateralization is known to differ between genders as well as age. Note that while a weak but systematic asymmetry between the hemispheres (difference between the strength of corresponding connections in the left and right hemispheres (Wilcoxon rank sum test, p < 0.01)) was observed in our data set in more than 20 percent of the connections —most prominently in the temporal cortex— we did not observe a significant variation of the asymmetry with gender, age, or handedness. This can probably be attributed to the relatively small size of our sample that provides information about the general structure and asymmetry of the brain, but not about its subject-specific variability.

Each node corresponding to a brain region is modeled by the FitzHugh-Nagumo (FHN) model, a paradigmatic model for neuronal spiking [30,31]. 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. Thus, the dynamics of the network reads

$$\begin{eqnarray} \epsilon \dot{u}_k =\,& u_k - \frac{u_k^3}{3} - v_k \nonumber\\ & {}+\sigma \sum_{j \in N_{\text{H}}} A_{kj} [ B_{uu}(u_j - u_k) + B_{uv}(v_j - v_k)]\nonumber\\ & + \varsigma \sum_{j \notin N_{\text{H}}}A_{kj} [ B_{uu}(u_j\! -\! u_k)\!+\! B_{uv}(v_j\!-\! v_k)],\qquad \end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray} \dot{v}_k =\,& v_k + a \nonumber\\ & {}+ \sigma \sum_{j \in N_{\text{H}}} A_{kj} [ B_{vu}(u_j - u_k) + B_{vv}(v_j - v_k)] \nonumber\\ & {}+ \varsigma \sum_{j \notin N_{\text{H}}} A_{kj} [B_{vu}(u_j\!-\! u_k)\! +\! B_{vv}(v_j\! -\! v_k)],\qquad \end{eqnarray} \tag{ 1b }$$

$k \in N_{\text{H}}$

$N_{\text{H}}$

$\epsilon = 0.05$

$(|a| > 1)$

$(|a| < 1)$

$a=0.5$

$$\begin{equation} \mathbf{B} =\left(\begin{array}{l@{\quad}c} B_{uu} & B_{uv} \\ B_{vu} & B_{vv} \end{array}\right)=\left(\begin{array}{l@{\quad}c} \cos\phi & \sin\phi \\ -\sin\phi & \cos\phi \end{array}\right)\!, \end{equation} \tag{ 2 }$$

with coupling phase $\phi = \frac{\pi}{2} - 0.1$, causing primarily an activator-inhibitor cross-coupling. This particular scheme was shown to be crucial for the occurrence of chimera states in ring topologies [32] as it reduces the stability of the completely synchronized state.

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_k complete rotations were realized. Throughout the paper we use $\Delta T = 5000$. Furthermore, we introduce hemispheric measures that characterize the degree of synchronization of the sub-networks and give complementary information. First, the hemispheric mean phase velocity is

$$\begin{equation} \langle \omega \rangle_H = \frac{1}{45}\sum_{k \in N_{\text{H}}}\omega_k, \end{equation} \tag{ 3 }$$

where H denotes either the left $(H=L)$ or the right $(H=R)$ hemisphere. Thus, $\langle \omega \rangle_H$ corresponds to the mean phase velocity averaged over the left or right hemisphere, respectively. To quantify the dynamical difference between the left and the right hemisphere we use the difference between these hemispheric mean phase velocities $\Delta \omega = \langle \omega \rangle_R - \langle \omega \rangle_L$.

Second, the hemispheric Kuramoto order parameter

$$\begin{equation} R_H(t) = \frac{1}{45} \left|\sum_{k \in N_\text{H}} \exp[i \theta_k(t)]\right| \end{equation} \tag{ 4 }$$

is calculated by means of an abstract dynamical phase $\theta_k$ that can be obtained from the standard geometric phase ${\tilde{\phi}_k(t) = \arctan(v_k/u_k)}$ by a transformation which yields a constant phase velocity $\dot{\theta}_k$. For an uncoupled FHN oscillator the function $t(\tilde{\phi}_k)$ is calculated numerically, assigning a value of time $0<t(\tilde{\phi}_k)<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(\tilde{\phi}_k)/T$, which yields $\dot{\theta}_k = \text{const}$. Thereby identical, uncoupled oscillators have a constant phase relation with respect to the dynamical phase. Fluctuations of the order parameter R_H caused by the FHN model's slow-fast time scales are suppressed and a change in R_H indeed reflects a change in the degree of synchronization. The Kuramoto order parameter may vary between 0 and 1, where $R_H=1$ corresponds to complete phase synchronization, and small values characterize spatially desynchronized states. Finally, we use the spatial correlation coefficient introduced by Kemeth et al. [33]:

$$\begin{equation} g_0(t) = \sqrt{\int_0^{\delta} g(t, D)\mathrm{d}D}, \end{equation} \tag{ 5 }$$

that measures the relative amount of synchronized oscillators. It is defined in terms of the normalized probability density function g(t, D), which is calculated as the probability of finding a distance D among all pairwise distances $\{D_{kj}\}= \{| e^{i\theta_j} - e^{i\theta_k}|\}$ between the states of all oscillators k, j, and generalizes the local curvature in systems with a spatial dimension. This distance is calculated using the dynamical phase on the unit cycle where the maximum distance of two oscillators is $D_{\textrm{max}} = 2$. For complete phase synchronization the distance between each pair of oscillators vanishes, i.e., $D = 0$ and $g(t, D) = \delta(D)$, hence $g_0(t)=1$, while a totally incoherent system gives a value of $g(t,0) = 0$, hence g₀(t) is small. Two oscillators are considered spatially correlated if their distance is smaller than some threshold $\delta = 0.01 D_{\max}$. The square root in eq. (5) arises because by taking all pairwise distances, the probability of oscillators k and j both being in the synchronous cluster is proportional to the square of the number of synchronous oscillators.

We investigate dynamical asymmetries emerging from the slight structural asymmetry of the brain hemispheres. Figure 2 shows how the different measures lead to the observation of a dynamical asymmetry with respect to the hemispheres of the average empirical structural brain network. Figure 2(a) displays the node-wise mean phase velocity $\omega_k$ for an intermediate coupling strength $\sigma = \varsigma = 0.7$ with random initial conditions. Note that the oscillators split into two visually well-distinguishable communities that coincide with the hemispheres of the brain network and have different hemispheric mean phase velocities $\langle \omega \rangle_H$. The left and right hemispheric mean phase velocities $\langle \omega \rangle_L$ and $\langle \omega \rangle_R$ and their difference $\Delta \omega$ vs. $\sigma = \varsigma$ are displayed in fig. 2(c) and (d), respectively. The values are calculated for one hundred different coupling strengths with $0 < \sigma \leq 1$ and step-size $\Delta \sigma = 0.01$. For every coupling strength an average over ten simulations with different sets of random initial conditions is plotted. For coupling strength $\sigma={\varsigma \ge 1}$ the system enters the frequency-synchronized regime, while phase synchronization measured by the Kuramoto order parameter sets in only later at $\sigma=\varsigma \approx 4.85$. Note that frequency synchronization means that all mean phase velocities $\omega_k$ are equal, but the phases may be different, while phase synchronization requires additionally that all phases are the same. It turns out that the difference $\Delta \omega$ assumes a maximum at $\sigma \approx 0.7$ and subsequently decreases again as both hemispheres enter the frequency-synchronized regime. However, these differences between left and right hemisphere do not imply different dynamical regimes in the sense of a partial synchronization pattern consisting of a desynchronized and a synchronized hemisphere, like in a chimera pattern. Nevertheless it can clearly be concluded that the network dynamics reflects the slight structural asymmetry. Figure 2(b) depicts the critical coupling strength for the transition between incoherence and frequency synchronization for a wider range of parameters in the ($\sigma,a$)-plane by a solid line for the coupled network with $\sigma=\varsigma$, and by a dashed line for the isolated hemispheres $(\varsigma = 0)$. The additional coupling between the hemispheres leads to a higher threshold value $\sigma_c$ for frequency synchronization.

Zoom In Zoom Out Reset image size

So far, we have used an averaged empirical matrix to detect a dynamical asymmetry. For a deeper insight it is important to consider all twenty available empirical structural brain networks individually. In all, we observe one of three transition scenarios from incoherence $(0 < \sigma < 1)$ to frequency synchronization ($\sigma > 1$ for $a=0.5$) with increasing coupling strength, as shown in fig. 3. They are distinguished by the difference of the hemispheric mean phase velocities $\Delta \omega$ exhibiting either a pronounced single maximum (d), or a (negative) minimum followed by a pronounced maximum (e), or essentially no dynamical asymmetry at all (f). However, in most cases (17 out of 20) a dynamical asymmetry was measurable by means of $\Delta \omega$.

Zoom In Zoom Out Reset image size

In the following we analyze to which extent the dynamical asymmetry can be attributed to the structural asymmetry of the network by introducing a structural asymmetry parameter ρ with $0 \le \rho \le 1$ that allows for a continuous tuning between the original structural brain network and a fully symmetrized network, in the sense that both hemispheres are identical. We introduce the coupling matrix elements of a network interpolating between asymmetric and symmetric hemispheres by

$$\begin{equation} A^*_{kj} = \rho A_{kj} + (1-\rho)\overline{A}_{kj}, \qquad \rho \in [0,1] \end{equation} \tag{ 6 }$$

with $\overline{A}_{kj} = \frac{1}{2} (A_{kj} + A_{k+45, j+45})$, where all indices are taken modulo 90. The resulting matrix $\{A^*_{kj} \}$ describes identical hemispheres if $\rho = 0$ and coincides with the original empirical matrix if $\rho = 1$. We observe that the dynamical asymmetry, expressed by the hemispheric difference of the mean phase velocities $\Delta \omega$, builds up as the structural asymmetry parameter increases (fig. 4(a)). The dynamical asymmetry is most pronounced for an intermediate coupling strength, not too small, but also not too close to the threshold of frequency synchronization. By integrating

$$\begin{equation} W = \int_0^1 \mathrm{d}\sigma \Delta \omega \end{equation} \tag{ 7 }$$

we obtain a dynamical asymmetry parameter W which is indeed almost linearly correlated with the structural asymmetry parameter ρ. The Pearson correlation coefficient of these two measures is $r_{\rho, W} = 0.96$. The dynamical asymmetry parameter increases linearly with ρ up to a certain degree of structural asymmetry, but then saturates and does not increase further. This means that even the slightest structural asymmetry results already in a slight dynamical asymmetry, i.e., there is no threshold behavior. However, a slight dynamical asymmetry here does not induce an immediate symmetry breaking as known from critical phenomena. The increase of dynamical asymmetry instead first increases linearly with the structural asymmetry. Beyond this regime of linear response, the dynamical asymmetry does not increase further if the structural asymmetry increases beyond a certain degree, and the real empirical structural asymmetry seems to be just closely above that value corresponding to saturation of sensitivity. These results are consistent with our hypothesis that both unihemispheric and bihemispheric sleep can be possible dynamical states of the same network.

Zoom In Zoom Out Reset image size

To achieve partial synchronization patterns we consider the inter-hemispheric coupling strength ς as an independent parameter that allows us to reduce the coupling between the hemispheres. This is motivated by the presumption that sleeping with one hemisphere at a time requires a certain degree of hemispheric separation [1]. All other parameters remain unchanged.

We analyze the parameter regime where the previously used average empirical structural brain network with identical inter- and intra-hemispheric coupling strength $\varsigma = \sigma$ exhibits qualitatively different behavior from that with separated hemispheres $\varsigma = 0$, i.e., these two cases correspond to different dynamical regimes. For both cases we have numerically determined the critical intra-hemispheric coupling strength $\sigma_c$ for which the system engages into the frequency-synchronized regimes, see fig. 2(b). As $\varsigma = 0$ leaves us with two disconnected sub-networks, these sub-networks are naturally easier to synchronize. Note that these two disconnected sub-networks technically result in two different critical coupling strengths. However, the difference between these critical values is very small and thus negligible. Consider a coupling strength σ that lies within the shaded area of fig. 2(b). There, a phase transition with increasing ς must be expected, since the system is frequency synchronized if $\varsigma = 0$, and completely incoherent if $\varsigma = \sigma$. We find that the frequency-synchronized solution indeed breaks down in one hemisphere. This gives rise to the partial synchronization pattern shown in fig. 5 where the left hemisphere is incoherent while the right is frequency synchronized, except for three small brain regions (hippocampus, gyrus parahippocampalis and amygdala). This shows up in the space-time plot, in the mean phase velocity profile, and in the hemispheric Kuramoto order parameter (although there is no perfect phase synchronization resulting in R_R < 1). Note that the incoherent, left hemisphere occasionally exhibits a high degree of synchronization that, in contrast to the right hemisphere, is unstable and vanishes after a short while. In general, partial synchronization patterns where different dynamical regimes occur in the two hemispheres can be found whenever a phase transition with respect to ς is expected, i.e., in the shaded region of fig. 2(b).

Zoom In Zoom Out Reset image size

Finally, we analyze the transition from frequency to phase synchronization which occurs at much higher coupling strengths than shown in fig. 2(b), e.g., for $a=0.5$ at $\sigma=\varsigma=4.8$ and at $\sigma=3.4$ for $\varsigma=0$. Here we use the temporal mean of the spatial correlation coefficient g₀ that is suitable to distinguish between phase $(g_0=1)$, frequency synchronization (g₀ < 1), and complete incoherence $(g_0=0)$. In contrast to the Kuramoto order parameter R, g₀ provides an arbitrary threshold δ and thereby gives a more pronounced transition from frequency to phase synchronization. Figure 6 shows that in the $(\sigma,\,\varsigma)$-plane of coupling strengths there exist regimes where a high degree of phase synchronization in one hemisphere coincides with a low degree of phase synchronization in the other (i.e., only frequency synchronization). Further, we find that the degree of phase synchronization expressed by g₀ may exhibit non-monotonic behavior as a function of ς . To a certain amount this is expected as we have seen before that coupling two hemispheres $(\varsigma \ne 0)$ decreases the degree of synchronization, we thus expect a maximum of g₀ at $\varsigma = 0$. However, in certain regimes of σ we find a subsequent maximum at $\varsigma \approx 2$ which implies that increasing the coupling between the hemispheres does not necessarily decrease the degree of phase synchronization. Furthermore, we find that this subsequent maximum is in principle possible in both hemispheres (indicating bistability, cf. panels (a) and (b)), but is realized by only one hemisphere at a time. A high degree of phase synchronization in one hemisphere thus suppresses phase synchronization in the other. This could be an important mechanism for the occurrence of unihemispheric sleep and should be investigated further. It is interesting that the two hemispheres can exchange their roles as one being phase synchronized and the other being only frequency synchronized, depending upon the initial conditions. For very strong σ and ς complete phase synchronization of both hemispheres is found (top right corners in fig. 6).

Zoom In Zoom Out Reset image size

We have investigated the dynamical asymmetry arising from the structural difference between the two brain hemispheres. It has been found that during the transition from complete incoherence to frequency synchronization an asymmetry regarding the degree of synchronization builds up, which can be quantified by the different mean phase velocities averaged over each hemisphere. We have shown that this asymmetry can be attributed to the structural asymmetry of the hemispheres, by introducing an asymmetry parameter which can interpolate between the empirical asymmetric brain network and an artificially symmetrized network. Furthermore, we have varied the inter-hemispheric coupling strength, while keeping the intra-hemispheric coupling strength fixed, to increase the degree of inter-hemispheric separation, ranging from isolated to fully coupled hemispheres. This has resulted in the observation of partial synchronization patterns similar to spontaneously synchrony-breaking chimera states. We have demonstrated that these partial synchronization patterns occur for coupling strengths where the isolated hemispheres are frequency synchronized while the brain network with equal intra- and inter-hemispheric coupling remains completely incoherent. By tuning the coupling between the hemispheres we have shown that at intermediate inter-hemispheric coupling one hemisphere becomes incoherent, giving rise to a chimera-like partial synchronization pattern.

These results are in accordance with the assumption that unihemispheric sleep requires a certain degree of inter-hemispheric separation. Moreover, it is known that the brain is operating in a critical state at the edge of different dynamical regimes [34], exhibiting hysteresis and avalanche phenomena as seen in critical phenomena and phase transitions [35–37]. By choosing appropriate coupling parameters, we have reported an intriguing dynamical behavior regarding the transition from frequency to phase synchronization. We observe that in this regime our brain model exhibits spontaneous symmetry breaking and bistability, where each hemisphere may engage into either of two dynamical states, characterized by a relatively high and low degree of synchronization. However, a high degree of synchronization in one of the hemispheres always coincides with a low degree of synchronization in the other. To sum up, the structural asymmetry in the brain allows for partial synchronization dynamics, which may be used to model unihemispheric sleep or explain the mechanism of the first-night effect in human sleep.

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer - 163436311 - SFB 910 and by project No. LO1611 with financial support from the MEYS under the NPU I program. LR is grateful to P. Jiruska and J. Hlinka for their hospitality and discussion during his stay in Prague. We are grateful to Antonín Škoch for preparing the example structural connectivity matrices.