Directed Flow of Information in Chimera States

We investigated interactions within chimera states in a phase oscillator network with two coupled subpopulations. To quantify interactions within and between these subpopulations, we estimated the corresponding (delayed) mutual information that -- in general -- quantifies the capacity or the maximum rate at which information can be transferred to recover a sender's information at the receiver with a vanishingly low error probability. After verifying their equivalence with estimates based on the continuous phase data, we determined the mutual information using the time points at which the individual phases passed through their respective Poincar\'{e} sections. This stroboscopic view on the dynamics may resemble, e.g., neural spike times, that are common observables in the study of neuronal information transfer. This discretization also increased processing speed significantly, rendering it particularly suitable for a fine-grained analysis of the effects of experimental and model parameters. In our model, the delayed mutual information within each subpopulation peaked at zero delay, whereas between the subpopulations it was always maximal at non-zero delay, irrespective of parameter choices. We observed that the delayed mutual information of the desynchronized subpopulation preceded the synchronized subpopulation. Put differently, the oscillators of the desynchronized subpopulation were 'driving' the ones in the synchronized subpopulation. These findings were also observed when estimating mutual information of the full phase trajectories. We can thus conclude that the delayed mutual information of discrete time points allows for inferring a functional directed flow of information between subpopulations of coupled phase oscillators.

Chimera states can be considered patterns of localized synchronization. As such they may contribute to the coding of information in a network. This is particularly interesting since systems with chimera states typically display multi-stability, i.e., different chimera configurations may co-exist for identical parameters [41][42][43][44]. Such a network may hence be able to encode different stimuli through different chimera states without the need for adjusting parameters or altering network structure. It is even possible to dynamically switch between different synchronization patterns, thus allowing for dynamic coding of information [45,46].
The mere existence of different levels of synchronization in the same network can also facilitate the transfer of information across a network, especially between subpopulations. Coherent oscillations between neurons or neural populations have been hypothesized to provide a communication mechanism while full neural synchronization is usually considered pathological [24,47,48]. A recent study reported chimera-like states in neuronal networks in humans, more specifically in electro-encephalographic patterns during epileptic seizures [49]. Admittedly, our understanding of these mechanism is still in its infancy also because the (interpretations of) experimental studies often lack mathematical rigor, both regarding the description of synchronization processes and the resulting implications for the network's capacity for information processing. What is the information transfer within and between synchronized and desynchronized populations? We investigated the communication channels between two subpopulations in a system of coupled phase oscillators. Since we designed that system to exhibit chimera states, we expected non-trivial interactions between the subpopulations. To tackle this, we employed the delayed version of mutual information which measures the rate at which information can be sent and recovered with vanishingly low probability of error [50]. Furthermore, assuming symmetry of our configuration, any directionality found in the system should be regarded of functional nature. With this we sought not only to answer the aforementioned question but to contribute to a more general understanding of how information can be transferred between subpopulations of oscillators. In view of potential applicability in other (experimental) studies, we also investigated whether we could gain sufficient insight into this communication by looking at the passage times through the oscillators' Poincaré sections rather than evaluating the corresponding continuous-time data. Such data may resemble, e.g., spike trains in neurophysiological assessments.
Our paper is structured as follows. First, we introduce our model in the study of chimera states [32,[51][52][53]. It consists of two subpopulations that either can be synchronized or desynchronized. We generalize this model by including distributed coupling strengths among the oscillatory units as well as additive Gaussian white noise. We briefly sketch the conditions under which chimera states can exist. Second, we outline the concept of delayed mutual information and detail how to estimate this for 'event-based' data where we define events via Poincaré sections of the phase oscillators' trajectories. With this tool at hand, we finally investigate the flow of information within and between the two subpopulations, and characterize its dependency on the essential model parameters.

II. MODEL
We build on a variant of the noisy Kuramoto-Sakaguchi model [54], generalized to M subpopulations of phase oscillators [55]. The phase φ j,µ (t) of oscillator j = 1, . . . , N in population µ = 1, . . . , M evolves according to where ω j,µ is the natural frequency of the oscillator j in subpopulation µ. Throughout our study ω j,µ is drawn from a zero-centered Gaussian distribution with variance σ 2 ω . Oscillator j in population µ and oscillator k in population ν interact sinusoidally with coupling strength C kj,µν and a phase lag α. The phase lag α varies the interaction function between more cosine (α → π/2) or more sine like-behavior (α → 0) and can be interpreted as a (small) transmission delay between units [34]. The additive noise term dW j,µ (t) represents mean-centered Gaussian noise with variance (strength) σ 2 W , i.e., dW j,µ (t)dW k,ν (t ) t = σ 2 W δ jk δ µν δ(t − t ). For the sake of legibility, we restrict this model to the case of M = 2 subpopulations of equal size N 1 = N 2 = N . In the absence of noise and in the case of identical oscillators, uniform phase lag and homogeneous coupling, i.e., for σ W = 0, σ ω = 0, α µν := α, and σ C = 0, respectively, the aforementioned studies establishes corresponding bifurcation diagrams [32,51,52].
To generalize the model towards more realistic, experimental setting, we included distributed coupling strengths via a non-uniform, (statistically) symmetric coupling between subpopulations. This can be cast in the following coupling matrix: with block matrices C ζ , C η ∈ R N ×N . The self-coupling within the subpopulation and the neighbor coupling between subpopulations are represented by C ζ and C η , respectively. Each block C ζ (or C η ) is composed of random numbers drawn from a normal distribution with mean ζ (or η) and variance σ 2 C . By this, we can tune the heterogeneity in the network. We would like to note that, in general, the chosen coupling is only symmetric for σ C > 0 in the sense of a statistical average, and that the expected coupling values ζ and η in each block only can be retrieved by the mean in the limit of large numbers. Although the symmetry between the two subpopulations 1 ↔ 2 is broken and only preserved in a statistical sense in the limit of large numbers, we verified that chimera states appeared in both configurations. That is, subpopulation 1 is synchronized and subpopulation 2 is desynchronized (SD) and subpopulation 2 is synchronized and subpopulation 1 is desynchronized (DS), for a particular choice of parameters using numerical simulations. Another consequence of distributed coupling strengths is that only a very few oscillators may have experienced negative (inhibitory) coupling, which can be neglected.
Following [51], we parametrize the relation between strengths by A = ζ − η with ζ + η = 1 and the phase lag α µ,ν by β µ,ν = π 2 − α µ,ν which here can be kept homogeneous for the entire network, i.e., β µ,ν := β. By this, we may recover the results reported in [32,51] in the limit σ C → 0, under the proviso of σ W = σ ω = 0. In brief, increasing A from 0 while fixing α (or β := π/2 − α) yields the following bifurcation scheme: a "stable chimera" is born through a saddle-node bifurcation and becomes unstable in a supercritical Hopf bifurcation with a stable limit cycle corresponding to a "breathing chimera", which eventually is destroyed in a homoclinic bifurcation. For all parameter values, the fully synchronized state exists and is a stable attractor; see also Fig. 4 in Appendix A. Subsequent studies demonstrated robustness of chimera states against non-uniformity of delays (α µν = α) [33], heterogeneity of oscillator frequencies (σ ω > 0) [56], and additive noise [57]. Although nonuniform phase lags lead to less degenerate dynamics and give room for more complex dynamics [33,53], we restrict ourselves to the case of uniform phase lags without compromising the essential phenomenology of chimera states.
The macroscopic behavior, i.e. the synchronization of the two populations, may be characterized by complexvalued order parameters, which are either defined on the population level, i.e., Z µ := N −1 N j=1 e iθµ,j , or globally, i.e., Z = (2N ) −1 M µ=1 N j=1 e iθµ,j . As common in the studies of coupled phase oscillators, the level of synchrony can be given by R µ := |Z µ |. Thus, R µ = 1 implies that oscillators in population µ are perfectly phase synchronized (S), while for R µ < 1 the oscillators are imperfectly synchronized or de-synchronized (D). For our two-subpopulation case, full synchrony (SS) hence occurs when R 1 1 and R 2 1, and chimera states are present if R 1 < 1 and R 2 1 or vice versa. The angular order parameter, Φ µ := arg Z µ keeps track of the average phase of the (sub)population. Fluctuations inherent to the model may affect the order parameter as illustrated in Fig. 1(b) (please refer to section III A for the numerical specifications). We therefore always considered averages over time, R µ t , when discussing the stability of a state. In fact, in our model chimera states remain stable for relatively large coupling heterogeneity, σ C > 0 presuming σ ω > 0 and σ W > 0, as is evidenced by numerical simulations; see also Appendix A, Fig. 4. The perfect synchronization manifold with R = 1 cannot be achieved; see Fig. 1(b). Further aspects of these noisy dynamics will be presented elsewhere.
Adding noise and heterogeneity to the system may alter its dynamics. In the present work we concentrated where · t and · j denote averages over time and oscillators, respectively. Parameters are: on parameter regions characterized by the occurrence of dynamic states that did resemble stable chimeras, i.e., where · t denotes the average over a duration of T = 9·10 6 time steps (after removing a transient of T = 10 5 time steps). Fig. 4 in Appendix A provides an overview and explanation of how parameter points A and β were selected.

A. Simulations
For the numerical implementation we employed a Euler-Maruyama scheme with ∆t = 10 −2 for N = 128 phase oscillators per subpopulation that were evolved for T = 10 6 [54]. We varied the coupling parameter A and the phase lag parameter β, while we fixed the width (standard deviation) of the natural frequency distribution to σ ω = 0.01 and of the coupling distribution to σ C = 0.01. The additive noise had variance σ 2 W = 0.001.

B. Mutual Information
Mutual information I(X; Y ), first introduced by Shannon [58], is meant to assess the dependence between two random variables X and Y . It measures the amount of information that X contains about Y . In terms of communication theory, Y can be regarded as output of a communication channel which depends probabilistically on its input, X. By construction, I(X; Y ) is non-negative; it is equal to zero if and only if X and Y are independent. Moreover, I(X; Y ) is symmetric, i.e., I(X; Y ) = I(Y ; X) implying that it is not directional. The mutual information may also be considered as a measure of the reduction in the uncertainty of X(Y) due to the knowledge of Y (X) or, in terms of communication, the rate at which information is being shared between the two [50]. The maximum rate at which one can send information over the channel and recover it at the output with a vanishingly low probability of error is called channel capacity, c = max p(x) I(X; Y ). In networks with oscillatory nodes, the random variables are the degree of synchronization of different subpopulations, and the rate of information shared will depend on the state of the system for each of the subpopulations.
Mutual information can be defined via the Kullback-Leibler divergence or in terms of entropies, where H(X) is the entropy of the random variable X, i.e., H(X) = − dx p X (x) log p X (x) with p X (x) being the corresponding probability density.

C. Delayed mutual information
For time-dependent random variables, one may generalize this definition to that of a delayed mutual information. This variant dates back to Fraser and Swinney [59], who used the delayed mutual information for estimating the embedding delay in chaotic systems -in that context the delayed auto-mutual information is often considered a 'generalization' of the auto-correlation function [60]. Applications range from the study of coupled map lattices [61] via spatiotemporal [62] and general dependencies between time series [63] to the detection of synchronization [64]. With the notation I XY := I(X; Y ) and H X := H(X), the delayed mutual information for a delay τ may read which has the symmetry I XY (τ ) = I Y X (−τ ) [65]. With this definition, one can measure the rate of information shared between X and Y as a function of time delay τ . In fact, we are not particularly interested in the specific value of the mutual information but rather focus here on the time delay at which the mutual information is maximal. Hence, we define τ max := arg max τ I XY (τ ). A positive (negative) value of τ max implies that Y shares more information with a delayed (advanced) X. This means there is an information flow from X(Y) to Y (X).

D. Delayed mutual information between subpopulations
a. Estimates using continuous-time data. When the time-dependent random variables are continuous time series u µ (t) and v ν (t) associated to populations µ and ν, the delayed mutual information can be estimated from Eq. (4) for X(t) = u µ (t) and Y (t) = u ν (t), We determined the probability densities using kernel density estimators with Epanetchnikov kernels with bandwidths given through a uniform maximum likelihood cross-validation search. For our parameter settings (254 oscillators and 10 6 samples), the resulting bandwidths ranged from about 0.10 rad to 0.18 rad. This software implementation is part of the KDE-toolbox; cf. [66] and [67] for alternative schemes. b. Estimates using event-based time data. For the aforementioned event signals in the subpopulations 1 and 2, i.e., discrete time points are defined as passing moments through the respective Poincaré sections. The probability densities to incorporate when estimating the mutual information are densities of events, or densities of times. We implemented the probability estimates as follows. Let S µ be a set of event times, i.e., µ,i stands for the time of the m-th event of oscillator i in subpopulation µ. Then, the probability density for an event to happen at time t in subpopulation µ is p µ (t) = P (t ∈ S µ ) and the probability of an event to happen at time t in subpopulation µ and time t + τ in subpopulation ν is The delayed mutual information can be given as We again computed the probability densities using kernel density estimators [66] but now involving Gaussian kernels. We also adjusted the bandwidth selection to a spherical, local maximum likelihood cross-validation due to the sparsity of the data and the resulting bandwidths ranged from about 25 to 35 time units. These results appeared robust when using the aforementioned uniform search; again we employed the KDE-toolbox.

IV. RESULTS
To determine the directionality of the information flow in the network we computed the time lagged mutual information within the subpopulations, I 11 (τ ) and I 22 (τ ) and between them, I 12 (τ ) = I 21 (−τ ). The first results for the event data outlined in Section III E are shown in Fig. 2(a). Since the recurrence of events mimicked the mean frequency of the phase oscillators, the mutual information turned out periodic. As expected, we observed (average) values of the mutual information that differed between I 11 (red), I 22 (blue), I 12 (black). This relates to the difference in entropy of the subpopulations, with the less synchronized one (µ = 2) being more disordered. The latter, hence, contained more entropy. However, since we were not interested in the explicit values of I µν (τ ), we could rescale the mutual information to maximum value which allowed for a comparative view when zooming in to the neighborhood of τ ≈ 0; see panel (c) of Fig. 2. The off-zero peak of I max := I 12 (τ max ) in τ = τ max clearly indicated a directed information flow, i.e. there is a functional directionality despite the structural symmetry of our model.
When comparing the estimates for the mutual information obtained from the event data with those from the continuous-time we found little to no difference in peak location. That is, the positions τ max of the maximum peaks, I max andĨ max , were nearly identical using either method as shown in panels (b) and (d) Fig. 2. Thus, to study the effects of varying model parameters on I µν (τ ), our subsequent analysis may solely rely on the event data approach.
Varying the values of A and β revealed a strong relation between the location of the peak of the delayed mutual information τ max and the relative phase arg (Z 1 /Z 2 ) between the two subpopulations; see Fig. 3(a). This convincingly shows that our approach to analyze the eventbased data reveals important information about the (otherwise continuous) dynamics of the subpopulations, here given by the phases (arguments) of the local order parameters. By contrast, the relative strength of local synchronization |Z 1 /Z 2 | had little to no effect on τ max ; see Fig. 3(b). These dependencies were inverted when looking at the value of mutual information, I max /I SS max . Consequently, the value of mutual information was affected by relative strength of local synchronization |Z 1 /Z 2 |after all the more a subpopulation is synchronized, the lower the corresponding entropy. However, effects were small and probably negligible when transferring to (selected) experimental data. More details about the normalization factor I SS max are discussed in Appendix B.

V. DISCUSSION
The delayed mutual information from event data, here, the passing times through the oscillators' Poincaré sections, agreed with that of the continuous data. This is an important finding as it shows that studying discrete event time series allows for inferring information-related properties of continuous dynamics. This offers the possibility to bring our approach to experimental studies, be that in neuroscience, where spike trains are a traditional measure of continuous neuronal activity, or economics, where stroboscopic assessments of stocks are common practice. Here, we used this approach to explore the information flow within and between subpopulations of a network of coupled phase oscillators. This information flow turned out to be directed.
Mutual information is the traditional measure of shared information between the two systems [58]. Our findings relied on the introduction of a time delay τ , which readily allows for identifying a direction of flow of information. This concept is much related to the widely used transfer entropy [68,69]. In fact, transfer entropy is the delayed conditional mutual information [70,71]. It therefore differs from our approach primarily by a nor- malization factor. Since here we were not interested in the absolute amount of information (flow), we could simplify assessments and focus on delayed mutual information -adding the conditional part jeopardizes the approximation of probabilities, i.e., estimating transfer entropy in a reliable manner typically requires more data than our delayed mutual information. A detailed discussion of these differences is beyond the scope of the current study. For our model we found that the delayed mutual information between the synchronized subpopulation and the less synchronized one peaked at a finite delay τ max = 0. Hence, there was a directed flow of information between the two subpopulations. Since we found that τ max < 0 for I 12 (τ ), the direction of this flow was from the less synchronized subpopulation to the fully synchronized one, i.e. 2 → 1 or D → S. In fact, the delay largely resembled the relative phase between the corresponding local order parameters, as shown in panel (a) of Fig. 3. We could not only readily identify the relative phase by encountering event data only, but our approach also allowed for attaching a meaningful interpretation in an information theoretic sense. This is promising since event data are, as stated repeatedly, conventional outcome parameters in many experimental settings. This is particularly true for studies in neuroscience, for which the quest on information processing is often central. There, networks are typically complex and modular. The complex collective dynamics switches at multiple scales, rendering neuronal networks especially exciting when it comes to information routing [72]. As of yet, our approach does not allow for unraveling dynamics information routing. This will require extending the (time-lagged) mutual information to a time-dependent version, e.g., by windowing the data under study. We plan to incorporate this in future studies.

VI. CONCLUSION
Estimating the delayed mutual information based on time points at which the individual phases passed through their respective Poincaré sections allows for identifying the information flow between subpopulations of networks. If the network displays chimera states, the information flow turns out to be directed. In our model of coupled phase oscillators, the flow of information was directed from the less synchronized subpopulation to the fully synchronized one since the first preceded the latter. Our approach is a first step to study information transfer between spike trains. It can be readily adopted to static well-defined modular networks and needs to be upgraded to a time dependent version to be applied to real, biological data.

FUNDING
This study received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement #642563 (COSMOS).

AUTHOR CONTRIBUTIONS
ND conducted the simulations and analysis and wrote the manuscript; AD set up the study, incorporated the kernel density estimation and wrote the manuscript, DB contributed to the first concepts of the study and wrote the manuscript; EAM set up the study, conducted simulations and wrote the manuscript.
For the study of mutual information, we restricted ourselves to dynamic states most similar to the 'stable chimera states' reported in [51]. Due to the various noise sources inherent to the model, the system dynamics undergoes fluctuations, as described in Section II. This poses additional challenges that we addressed by examining specific time averages of the complex-valued order parameter in each subpopulation. Since fluctuations may cause the system to temporarily drift off the synchronized manifold, we may not simply consider |Z 1 | = 1, |Z 2 | < 1 to identify chimera states. Instead we require the order to be sufficiently different between the two subpopulations, i.e., we threshold |Z 2 | t / |Z 1 | t , defining the region R 1 shown in Fig. 4(a). Stability diagram based on analyzing complexvalued order parameters. Regions R1 (a) and R2 (b) in green outline parameter values (β, A) for which |Z2| t/ |Z1| t < 0.8 and max(|Z2|) t − min(|Z2|) t < 0.1, respectively. The states that are most similar to 'stable chimera states' [51] were identified within the intersection R1 ∩ R2 (outlined in red). The outlined regions are bounded by parameter regions known to exhibit stable chimeras in the limit of σω = σC = σ 2 W = 0 [51]; saddle node and Hopf bifurcation curves are shown as black solid and dashed curves, respectively. The other parameters have been fixed at σω = 0.01, σC = 0.01, and σ 2 W = 0.001.
Furthermore, we presume amplitude variations to be small. To that end, we thresholded max(|Z 2 |) t − min(|Z 2 |) t , defining the region R 2 shown in Fig. 4(b).
Here, · t denotes the average defined over the maxima computed over time windows of length T = 10 3 time steps, over a duration T = 9 · 10 5 after removing a transient of T = 10 5 time steps. Finally, chimera states with parameter values (β, A) ∈ R 1 ∩ R 2 were selected for further analysis of mutual information.
Appendix B: Normalization factor I SS max Our model displays a baseline dependency of the maximum peak of mutual information I SS max on the parameters β and A, even when the system is in the fully synchronized state SS, see Fig. 5. We therefore studied the dependency of the maximum peak of the delayed mutual information for chimera states normalized with I SS max , I max /I SS max . For every set of parameters we computed the delayed mutual information on two different solution that have been forced by changing the initial conditions; the maxima of each of these solutions (SS and SD) are I SS max := I SS max (A, β) and I max := I SD max (A, β). These are the quantities displayed in panels (c) and (d) of Fig. 3. The normalization factor given by the peak value of mutual information in a fully synchronized state SS, I SS max := I11(τmax) = I22(τmax) as a function of (β, A). The other parameters are kept fixed at σω = 0.01, σC = 0.01, and σ 2 W = 0.001.