Fractal basins as a mechanism for the nimble brain

An interesting feature of the brain is its ability to respond to disparate sensory signals from the environment in unique ways depending on the environmental context or current brain state. In dynamical systems, this is an example of multi-stability, the ability to switch between multiple stable states corresponding to specific patterns of brain activity/connectivity. In this article, we describe chimera states, which are patterns consisting of mixed synchrony and incoherence, in a brain-inspired dynamical systems model composed of a network with weak individual interactions and chaotic/periodic local dynamics. We illustrate the mechanism using synthetic time series interacting on a realistic anatomical brain network derived from human diffusion tensor imaging. We introduce the so-called vector pattern state (VPS) as an efficient way of identifying chimera states and mapping basin structures. Clustering similar VPSs for different initial conditions, we show that coexisting attractors of such states reveal intricately “mingled” fractal basin boundaries that are immediately reachable. This could explain the nimble brain’s ability to rapidly switch patterns between coexisting attractors.

Here we expand concepts with further details as they appear in the body of this work.Specifically, network concepts for complex systems, dynamical systems and fractals concepts, and machine learning concepts for neuroinspired models of the brain.

Networks and Complex Systems 1.Networks
Networks play an important role in many complex systems, including the brain.A networked representation of the human brain, where different regions are treated as vertices and interactions between them as edges, allows us to explore both the structural and dynamical properties.
In this work, we use the phrase network synonymously with the mathematical conception of a graph G = (V (G), E(G)) which is a set of vertices (nodes) V = {v 1 , v 2 , ..., v N }, with cardinality |V (G)| = N and edges (links) E ⊆ V × V with |E(G)| = M .The adjacency matrix A is a matrix representation of the graph with entries where w i,j is a weight.A simple graph is an unweighted graph with no selfloops or multiple edges.Thus, each entry [A] i,j ∈ {0, 1} and [A] i,j = 0 if i = j.

Subgraphs
Informally, a subgraph is the graph that results from erasure of some of the vertices and edges of the original graph.That is, a graph H = (V (H), E(H)) is called a subgraph of a graph G if V (H) ⊆ V (G) and E(H) ⊆ E(G) and |V (H)| = N ′ ≤ N , |E(H)| = M ′ ≤ M .Now suppose the vertices in G are relabeled so that the corresponding vertices in H have the same label, that is V (G) = {v 1 , v 2 , ...v N }, and V (H) = {v 1 , v 2 , ..., v N ′ }.Then we call H an induced subgraph of G if E(H) = (V (H) × V (H)) ∩ E(G).In other words, if for every vertex in H, the edges formed between each pair of those vertices in the graph G are also found in H, then H is an induced subgraph.
A graph is called a proper subgraph if it is a subgraph of G and it is not identically the same as G, i.e., it has strictly fewer nodes than G. Suppose a graph is partitioned into two proper induced subgraphs H, K ⊂ G such that V (H) ∪ V (K) = V (G) and V (H) ∩ V (K) = ∅; that is, H and K are mutually exclusive induced subgraphs.Then G is called connected if for every such H and K pair in G, E(H) ∪ E(K) ̸ = E(G).That is, if for every pair of proper induced subgraphs satisfying the properties above, there is at least one edge between the pair, then G is connected.Said another way, there is a path from every vertex to every other vertex.If G is not connected, then it is called disconnected.

Graph Symmetries
Symmetry is an important property of graphs that feature prominently in the topic of dynamics on networks, specifically synchrony, though it is not a necessary condition [42].By symmetry, we mean an automorphism of a graph.A graph automorphism is a permutation of the vertices which leaves the adjacencies unchanged, meaning is "essentially" the same graph.Stated in terms of the adjacency matrix A of a graph, there is a similarity transformation by a permutation matrix P such that A = P T AP , where P can be constructed from the identity matrix by permuting rows and columns only, also here P T is the transpose of P .Such an automorphism is called a symmetry of a graph if it exists and the graph is called symmetric.Note that the term symmetric here does not mean a symmetric adjacency matrix, but the phrase is usually only used when there exists such a non-trivial permutation.Also, P ̸ = I, the identity , which otherwise always works.

Diffusion Tensor Imaging (DTI) Graphs
Diffusion Tensor Imaging (DTI) is a magnetic resonance imaging technique used to find the structural organization of white matter tracts in the human brain.Since its formation in 1994, the technique has been widely used in brain research.The DTI data can be processed to estimate the brain fiber tracts using a technique called tractography.For the graph shown in Fig. 1, which is Subject 1 from [9], a process called deterministic tractography, which is available in a software package [45] called the diffusion toolkit, was performed.The brain was parcellated into 83 regions of interest (ROIs) and fiber tracts from each were computed.Using the ROI's as individual vertices and fiber tracts between ROIs as edges, a weighted structural brain network was generated.The DTI network is then converted to a simple graph as seen Figure 1: Network is formulated from diffusion tensor imaging (DTI), which is a technique for inferring white matter fiber tracts in the brain.In this case there are 83 regions of interest (ROI), and thus 83 vertices in the corresponding network structure.The size of each node scales to their betweenness centrality.This network has no non-trivial symmetries.
in Fig. 1 by changing all non-zero weights to 1 as was done in [13].When we refer to a DTI network throughout the main text and in the SI, it is this graph that we are referring to.
2 Dynamics On Networks, Synchronization and Chimera States

Synchronization
In a coupled dynamical system, synchronization means a 'rhythm' or 'a kind of sympathy' between two or more of its dynamical entities.The simultaneous firing of neurons in the human brain, synchronized flashing of light by fireflies, and synchronized clapping by an audience are some of the popular examples [34].Synchronization arises due to interactions between dynamical elements: while the individual units try to move according to their own set of rules, attractive interactions or coupling with others tries to bring harmony between them, which in turn leads to synchronization.Below, different forms of synchronization are discussed, each of which serve a role in our analysis of patterns in the brain, as deduced by our Vector Pattern States (VPS) to distinguish basins.

Identical Synchronization
The most basic and perhaps standard version of synchronization is called "identical synchronization" (IS) or "complete synchronization" (CS).In a complex system of coupled dynamically evolving units, this state occurs when all units asymptotically "beat together", identically.
Coupled elements are each modeled as differential equations with a vector field defined by a function f described on a graph as vertices, with the edges describing communication between the vertices, through a coupling function h.A generic model used for study in such systems in continuous time is written as: Here, ẋi represents the derivative of x i ∈ R d with respect to time, the pair ẋi = f i (x i ) represents the isolated dynamics associated with an individual vertex, [A] i,j are entries of the adjacency matrix and h : R d → R d is the coupling function between a pair of vertices.In the case where we have f i = f (∀i) (meaning that the uncoupled oscillators obey identical dynamics) and h(x i , x j ) = 0, when x i = x j a synchronization manifold can be defined as Any collective state of all units which occurs in M is called identically synchronous.Generally, states which asymptotically converge to the synchronization manifold are of interest, initial conditions which do so are defined as being in the basin of synchronization.Thus the set, of all initial conditions which are asymptotically attracted to the synchronization manifold is known as the synchronization basin.

Stability of identical synchronization
In their seminal work relating to identical synchronization, Pecora and Carroll [33] analyzed the (local) stability of the synchronous state of a network coupled system.Let ⃗ x(t) = (x 1 (t), x 2 (t), ..., x n (t)), then a synchronous state is called stable if starting with ⃗ x(t) ∈ M, there exists a δ > 0 and That is the synchronous state is called stable if an infinitesmal perturbation to that state away from the synchronization manifold returns to M asymptotically.In [33] details of the conditions under which an identical system will have a stable synchronous state are stated.

Phase Synchronization
Phase synchronization between two signals x i (t) and x j (t) is a generalization of identical synchronization.Whereas for identical synchronization, x i (t) → x j (t), as t → ∞, for phase synchronization, a slightly weaker but otherwise similar condition must hold after one of the signals is shifted by a phase.If there exists a phase shift, τ > 0 such that ∥x i (t) − x j (t − τ )∥ → 0 as t → ∞, then the two signals are defined as phase synchronous.This property is prominent in our analysis of brain patterns.

Generalized Synchronization
Gneralized synchronization is another weakening of the concept of synchrony.Two signals x i (t) and x j (t) are defined as generalized synchronous if ∥x i (t) − f (x j (t))∥ → 0 for some nonidentity function f (x) ̸ = x.In the identity case that f (x) = x, then this would simply be a statement of identical synchronization.Usually the statement is in terms of a smooth function f , and the set of points D = (x, f (x)) on which orbits (x i (t), x j (t)) are attracted is called the synchronization manifold, and this phrase can be used even in the identical case.It is certainly possible for a system to exhibit phase and generalized synchrony.

Approximate Synchronization
If two signals x i (t) and x j (t) approach a synchronous state, such as identical synchrony, but they do not converge to that state, then the synchronous state is stable but not asymptotically stable.In other words, if for some time T > 0, and small ϵ > 0, then ∥x i (t) − x j (t)∥ < ϵ for all t > T , but they do not converge to zero, then the signals are approximately (identically) synchronous.

Chimera States
A chimera state is defined by the presence of both synchronous and asynchronous behaviors in a complex system.Specifically, if a collection of coupled dynamical systems with corresponding signals, {x 1 (t), x 2 (t), ..., x n (t)} results in at least one subset of these that synchronizes (in terms of one of the types of synchronization noted above), and a distinct other subgroup that does not synchronize, then this is called a chimera state.Amongst those that synchronize, there may be just one large synchronous cluster, or perhaps several clusters that synchronize but not with each other, this is still called chimera, as long as there is also a subset that is asynchronous.13) in each of 8 clusters.The x-axis represents vertex pairs while the y-axis contains all initial conditions forming forming a cluster.The color intensity represents entries of e l .(i) The color-coded network shows the chimera state in (f).The vertices with the highest degree of synchronization are colored white while the remaining ones are shown in black.(j) shows x i versus time for this chimera state and is plotted as follows: first, x i (t) and x j (t − τ * i,j ) for nodes i and j having the smallest L(i, j, τ * i,j ) value is plotted; thereafter, the phase-shift for any new node k is decided such that L(i, k, τ * i,k ) minimum over all the pre-existing nodes i.
Chimera states have long fascinated the synchronization community, be-ing first named chimera states (in an identical synchronization context) by Abrams and Strogatz in 2004 [2] though this phenomenon was reported earlier by Kuramoto and Battogtokh [20].Though originally the term chimera state was reserved for identical oscillators, the literature is now filled with examples in which the term chimera state is used for non-identical oscillators [21,38,27,18,22] and so we have chosen to use the term in the broadest sense.The definition used here allows states that are sometimes called cluster synchronization or simply a clustered state [19,28,46,6,7], however it has been observed that a change in the coupling strength is enough to create a chimera state by isolated desynchronization of all clusters except for one [32].An example of one such representation of a chimera state can be found in Fig. 2(j).
The kind of synchrony present defines the details of the chimera state, whether the synchronous group be identical, generalized, phase synchronous, or approximate.In our synthetic model of HR oscillators on the DTI network as well as smaller models, we observe all of these.The stability of synchrony associated with a chimera state has some dependence on symmetry [32].

Nearly Synchronous States
Expanding upon the concept of approximate synchrony, Sec.2.5, here for the HR model, we synonymously state as nearly synchronous.Consider Fig. 3 (a) and (b) in the model with only electrical coupling, Eq. 2 with dynamics f given by Eq. 19 and coupling function h from Eq. 20, yield identical synchronization as shown in Fig. 3 (a).However, in Fig. 3 (c) we have added a chemical coupling term changing the coupling function h to the one from Eq. 21, which is the model used by [18] (see Sec. 4).As can be seen, while the neurons still appear to be synchronized in phase, the peaks no longer match.In Fig. 3 (e) we have included both the electrical and chemical coupling terms and additionally we have introduced a slight parameter mismatch.The model parameter a, which was previously set to 1 for all oscillators is now drawn from a normal distribution a ∼ N (1, 0.1).This is an example of nearly synchronous states of which is possible to examine the linear stability [44], even in the context of chimera [39].The HR model with the same coupling parameters as (c) but with the model parameter a drawn from a normal distribution a ∼ N (1, 0.1).(f) Zoom in panel of (e) for better visualization.In this case, the oscillators peak in the same general time frame but these spikes have different peaks and different shapes.This represents a more generalized synchronization pattern.
3 Basin Structure and Fractals

Fractal Dimension
Examples of fractals are abundant in nature: they exhibit patterns across many scales that are self-similar, in some form.Real world examples, where the self-similarity is approximate and fine details do not continue indefinitely, include the coastlines of Britain (fractal dimension ≈ 1.21 ) and Norway (fractal dimension ≈ 1.52) [12].Newton's fractals, Von-Koch curves, Cantor sets, Sierpinski's triangle, and the Koch Snowflake [12] are examples of mathematical fractals that continue indefinitely.Fractals can be understood as structures with roughness that can be represented by a real number, called their dimension.A definitive property of fractals is the concept of non-integer dimensionality [12], as this fractional dimension motivates the very naming of sets as "fractals".Whereas fractals are all defined on scaling of data density relative to an exponent describing the dimensionality, manifold dimensionality describes the integer number of coordinates necessary uniquely parameterize a point.For instance, a straight line is 1 dimensional, a plane surface is 2 dimensional, and a cube is 3 dimensional.Fractal dimensionality is generally a description of the scaling behavior of a set with respect to decreasing window perspective; a non-integer fractal.It is a concept that applies to smooth as well as "rough" sets.From this perspective, the coastline of Norway is rougher than Britain, and correspondingly, its dimension is also farther from nearby integer numbers 1 and 2 [23,5,4,40].
The fractal dimension is formally defined by the Hausdorff dimension, and various popular methods serve as useful estimates.For example, the similarity dimension, the correlation dimension, and the box dimension are each useful data driven estimators of the fractal dimension with certain underlying assumptions on the data.For practical computational reasons, here we use the box dimension method to calculate the fractal dimension [35,15,37].

Box Dimension
The box counting dimension [15] is a popular estimator for fractal dimension.Squares (boxes or hyperboxes in higher dimensions) of side length ϵ are chosen, and the number of 'boxes' (N (ϵ)) of this size needed to cover the entire set are counted.The box dimension (d) is understood as the scaling exponent, in the limit as the area (or volume) of the box goes to zero with: supposing the limit exists [43].If the above limit does not exist, one may still find the upper and lower box dimensions of the set.Efficient algorithms for estimating the box dimension exist [37].Fig. 4(a) shows an example of a boundary set obtained using Eq.20, while Fig. 4(b) shows the plot of ln( 1 ϵ ) vs ln(N ) for this boundary set.The slope of the line of best fit is interpreted as the limit given in Eq. 5, which in this example is equal to 1.27.As expected, the dimension of the boundary set is non-integer.Note that in practice due to finite data set size, as ϵ → 0, N (ϵ) saturates after a certain ϵ value is reached, which is due to the fact that N (ϵ) can not increase further once each point in the set is covered by a unique square.

Fractal Boundaries and Riddled Basins
Complexity can be found even in the structure of the basins of attraction in multi-stable systems.The feature of interest to this work is the possibility of the concept of a high level of "non-smoothness" of a fractal basin boundary when moving between sets of initial conditions attracted to one basin and those leading to a different one.And in particular, we described a scenario that arises for the HR system and some parameters where a particularly highly intermingled fractal basin boundary called riddled-basins occurs.We shall describe these scenarios.
In general, a boundary point of a set has the property that every neighborhood of the point, no matter how small, includes some points both in the set, and not in the set.The boundary set of a given set is thus the union of all such boundary points.The boundary of an open disc in the plane, for example, is a circle, which clearly is a smooth set.However, in dynamical systems it is not uncommon that boundary sets of basins of attraction can be fractals.The set of points forming the boundary between the basins of different attractors need not form a smooth curve, and instead this boundary set may have a fractal dimension [26].Our proposal herein is that the key feature of complex networked models of the brain with the observed behavior of rapid switching between various patterns is the presence of fractal basin boundaries [14,26].In particular, there is an apparently rich "intermingling" of these boundaries as the present phenomenon of what is called riddled basins [3,29,10,30]; see Fig. 7. Riddled basins occur in many physical settings, notably in the basin structure, for example, of magnetic potential wells [47], but the fact that we point this out in neurophysiology as shown in Fig. 7, as a crucial mechanism behind the ability to quickly switch focus, the nimble brain theory, is a unique and new interpretation of this now classical concept from dynamical systems.

Vector Pattern States
Here, we describe in more detail what we have introduced as Vector Pattern States (VPS) that allow comparisons between different initial conditions The purpose of the VPS is to classify the various fully synchronous, chimera, and incoherent states which may coexist in the basin.Thus, we associate which components of the complex network are in some form of synchrony, including the weaker approximate synchrony and phase shifted synchrony as presented in Section 2. This is a crucial technological step for mapping the basins as described in the next section.
In general, synchrony is defined asymptotically (i.e. as t → ∞), which is impractical.Instead for a finite time interval time series a similarity measure must be introduced to determine if a system is 'close enough' to be considered synchronized.Let I ⊂ R and ∥u∥ L 2 (I,M ) be the norm given by where M is a compact subset of R d and ∥ • ∥ is the Euclidean norm ∥u∥ 2 on R d .We abuse notation and also denote ∥u∥ L 2 (I,M ) by ∥u∥ 2 2 .Also, we denote u τ (t) = u(t − τ ) as the function u evaluated at t − τ .Other L p norms could be used to define the similarity measure, due to equivalence of L p norms in finite dimensional spaces.We will investigate the influence of different norms in future research.
An example of the similarity measure between the trajectories x i (t) and x j (t) of a pair (i, j) could be an average over This is fine in the CS case, however as we are interested in a more general synchronization case, the similarity measure should also allow for the possibility of a non-zero time shift τ between the time series.Thus, we consider the following similarity measure The time shift that minimizes Eq.9, τ * i,j , is of interest in this work since synchronization is viewed from a less restrictive criterion than CS, allowing for time shifts.For real-valued time series (d = 1), note that L(i, j, τ ) is related to cross-correlation, in the sense that is maximized at a particular τ , which happens to be the same τ * i,j that minimizes L(i, j, τ ).In other words, Thus, it is possible to find τ using the well-studied method cross correlation.In practice, not only is the time interval finite, the sampling rate is finite as well, allowing the integrals to be replaced by sums and the cross correlation can be found with with τ * i,j = arg max τ R x i ,x j (τ ) otherwise we say that τ * i,j = 0, and when In our simulations, the time interval and sampling rate are chosen so that Eq.( 12) gives a good approximation of τ * i,j .For instance, for the HR system we select the fast-variable (x-coordinate), so the time series has sufficiently large number of cycles for the chosen time interval and step size.
For the case of phase synchronization (where only phase difference is relevant), the τ * i,j for every i ̸ = j encodes all information about different final states, that is, they distinguish between chimera, complete synchronous, and completely incoherent states.In the generalized version of synchronization discussed here, not only are the τ * i,j important, but so is the difference between the 'profile' of the time series, so L(i, j, τ * i,j ) should be taken into consideration as well.Therefore the state of the pair (i, j) can be further characterized by the ordered pair (τ * i,j , L(i, j, τ * i,j )).We are now ready to define the vector pattern state (VPS) for all (i ̸ = j) e l = (τ * ,l 1,2 , τ * ,l 1,3 , . . ., τ * ,l N −1,N , βL l (1, 2, τ * ,l 1,2 ), βL l (1, 3, τ * ,l 1,3 ), . . ., βL l (N −1, N, τ * ,l N −1,N ), ( 13) where β > 0 allows for scaling the relative importance of the phase shift vs. the similarity measure between the time series.Each initial condition will lead to its own VPS.Though such states may not be unique, for instance, an initial condition which leads to complete synchronization will have the VPS: e l = (0, 0, 0, ..., 0) (in the limit as t → ∞), and thus each initial condition which leads to CS will have the same VPS.In the above mentioned case of phase synchronization, the VPS will be given by e l = (τ * 1,2 , τ * 1,3 , . . ., τ * N −1,N , 0, 0, . . ., 0) (in the limit as t → ∞).

Mapping the Basin Structure
In the last section we described our VPS so as to compare components of network coupled dynamics.With this, we can map the entire space of initial conditions to decide a basin structure.The VPS allows for characterization of states which may have both a phase shift as well as differences in the 'shape' of the time series.As discussed above, in the CS scenario, all components of the VPS will be ≈ 0. Similarly in the phase synchronization, all components corresponding to L(i, j, τ * i,j ) will be 0, while the τ * i,j themselves become bounded from above.However the utility of the VPS is best expressed in chimera states and generalized synchronization.When some of the vertices become asynchronous and others follow a CS state, the VPS will have some components nearly 0 and others large.
We examine the CS case in more detail, noting it is analogous to the more generalized synchronization.Finding the basin structure is straightforward if the pairwise distances (under whatever measure is chosen) between all oscillators are small enough.The system may be considered in CS and the initial condition leading to this state may be marked as one color.Any initial condition not leading to CS will then be assigned a different color.However, this becomes much more challenging in chimera states, since not all pairs will reach CS.It may seem that classifying all states which have at least one synchronized pair as well as at least one desynchronized pair as in the same basin will be sufficient, but chimera states range from being nearly incoherent to being nearly coherent and classifying all as the same may hide very rich structure.
To resolve this, each VPS (e l ) is stacked into a matrix, This matrix encodes all final states evolving from each initial condition, which are indexed by l.E is then clustered, in this case using k-means, with each cluster assigned a color in the basin plot.This leaves the choice of k, which was done using a variation on the elbow method.An example of how this is done is shown in Fig. 6.
In the elbow method after a cluster has been chosen from k-means, the value of W (C) from Eq. 17 becomes the dependent variable, with k as the independent variable.When examining this relationship across multiple k's generally there will be a discontinuity which is clearly visible, which will be the optimal value, though there are situations where this can fail [16].In a slight variation, to identify the elbow, a log-log plot is used instead, and the last k before the majority of values fall on the linear trend is selected as the 'elbow' as shown in Fig. 6 (b).
The effects varying k has on the basin structure are shown in Fig. 7, the initial conditions chosen from the random plane shown here were drawn from a 1500 × 1500 grid, sampled over the unit square ([0, 1] 2 ).The initial conditions for all other oscillators were initially chosen uniformly at random from the uniform distribution U(−1, 1) and were then held constant with the exception of the initial conditions drawn from the random plane (which were The complexity of the basin can be seen even for small values of k, the optimal k balances between too little and too much detail in the basin structure.In Fig. 7, k = 8 was found to be optimal using the version of the elbow method discussed above.Note that the basin used in this example differs from the one shown in Fig. 6, hence the different values of k.

Clustering and k-Means
Often insights into a problem can be drawn by clustering data together to infer meaningful relationships.Numerous methods exist to perform clustering, however k-means is a popular, robust and relatively computationally inexpensive method for clustering a 'cloud' of data points together.With data x ∈ R S×d , where S is the sample size and d is the dimension, and making the assumption that data is all of the quantitative type so that distance has meaning, then k-means is performed by minimizing the within cluster mean squared Euclidean distance between points.Choosing 1 ≤ k ≤ S, let C be a clustering C = {C 1 , C 2 , ..., C k }, with C k ⊂ {1, 2, ..., S}.Then we define the dissimilarity measure d to be the squared Euclidean distance [16]: and the within cluster mean where I(i ∈ C k ) and I is the indicator function.The within cluster squared error (W ) is defined as, Note that d(x i , xi,j ) is the distance between a data point and the centroid of the cluster it is in.The k-means problem is to minimize this within cluster Figure 7: Effect of varying k on the basin plots.Even for small k, a complex basin structure is observed.The choice of k can be seen as a trade off between too much and too little detail in this context.
squared error, that is to find the cluster Ĉ with As noted above, there are numerous clustering techniques, and one can even generalize k-means by designing the distance function d(x, y) to be a different dissimilarity measure than the more common squared Euclidean distance, such a clustering method is called k-medoids [16].

Networks of coupled oscillators
Here, we show that the construction of basins of attractions via VPS can also be applied to networks of coupled oscillators.Initially, we illustrate our approach for two small size populations that are globally coupled to each other [1].Then, we apply the VPS construction to a network that does not have full permutation symmetry.In both scenarios our approach is able to capture the basin structure, illustrating its applicability.

Coupled HR neurons
Hindmarsh and Rose (HR) developed a model of neuronal firing in the brain [17].We have used this model as a complex network of coupled oscillators as a semi-synthetic system using a true DTI network.Each HR neuron follows the following dynamics: where x = [x, y, z] T .Additionally two different coupling functions were used, the first one with, and the other one with For the standard coupling scheme with h 1 given in Eq. 20, the values of a, b, c, d, x R , and s are selected based on earlier research works, see [11,17].In this work, we set these parameters to a = 1, b = 3, c = 1, d = 5, x R = −0.5(1+ √ 5), s = 4, and I = 3.27.This will be referred to as the standard HR model.The value of x R is chosen from the condition ẋi = ẏi = 0 (i.e. the fast variables) while z i (the slow variable) is neglected.With this choice of parameters an isolated HR neuron exhibits a periodic time evolution.To decide if the time evolution is periodic or not, we follow the approach of Poincare [43].Fig. 5 highlights how the behavior of the system changes when the parameter r is varied.For a given set of parameters, the number of cuts a moving HR neuron makes on the Poincare plane are counted.The Poincare plane is any plane transverse to a time-continuous flow, which in this case is the y − z plane located at x = 0. Fixing I and decreasing r towards 0, the y coordinates of all intersections on the Poincare plane are shown, such that x(t) < 0 and x(t+dt) > 0, where dt denotes the integration time step.The resulting plot (Fig. 5) is called a bifurcation diagram and the type of bifurcation is a period-doubling bifurcation of limit cycles.A finite number of cuts (let's say N c ) indicate a periodic trajectory with period N c , while the limit N c → ∞ shows a chaotic evolution.
An alternative model, which includes chemical coupling in addition to the electrical coupling, that is with h 2 from Eq. 21, the parameters are: a = 1, b = 3, c = 1, d = 5, s = 4, r = 0.005, x R = −1.6,I = 3.25, σ = 0.5, α = 0.03, which is known to contain chimeras [18].Here the chemical and electrical adjacency matrices are assumed to be the same, unlike in [18].This will be referred to as the chemical coupling model.An example of the dynamics is shown in Fig. 3 (c) and (d), and a version of this model where the parameter a for each oscillator is drawn from a normal distribution is shown in Fig. 3 (e) and (f).We consider two coupled populations as illustrated in Figure 8 a)

Two populations of globally coupled oscillators
where θ k i is the phase of the i-th oscillator in population k ∈ {1, 2} and ω is the oscillator frequency.Based on earlier works [1,24,31], we suppose that σ 11 = σ 22 = µ > 0, and σ 12 = σ 21 = ν > 0, with µ > ν.Also, by rescaling time, we may set µ + ν = 1 and introduce useful parameters A = µ − ν and β = π/2 − α that determine the existence of chimeras states in the system.For identical oscillators, we can change coordinates to the rotating frame ω, θ k i → θ k i − ωt, so the natural frequency can be fixed at ω = 0.The system (22) has full permutation symmetry, allowing a low-dimensional description in terms of macroscopic variables in the thermodynamic limit (N k → ∞) [1].In particular, the authors in [24] use the resulting reduced system to describe the basins of attraction of chimera states.The basin structure corresponds to the infinite size system, but it is expected to remain roughly the same for large populations sizes [31].We tested our VPS construction (results not shown), reproducing a basin structure that contains a few distinct states, as predicted by the infinite size system [24].Here, we consider the interesting case when the population size is small, and the dynamics of the full system deviates from the evolution of the reduced system [31].In fact, in this scenario, the evolution still has a low-dimensional description but the macroscopic variables do not evolve in a closed form equation, so the reduced system becomes intractable analytically, see [25] for details.
We consider the case of N = 5, A = 0.1 and β = 0.025, where it is known that stable chimeras exist [1,24,31].For initial conditions we consider equally spaced phases over the circle, i.e., for a fixed ϑ ∈ [0, 2π] so the system starts close to a incoherence state.We consider ϑ = 0.1, transient time of 2000 time steps, and integrate over 1000 more time steps to obtain the final state, sampling every 0.02 time steps.The basin plot is performed in the plane (θ 1 1 (0), θ 2 1 (0)) for which the initial conditions are varied over the interval [0, 2π], see Figure 8 (b).We see that the basin has coexistence of different chimera states and (complete) phase synchronization state, as described in [31], and our method is capable of identifying these cases.13) and a grid with 1248×1248 uniformly sampled initial conditions.

Coupled oscillators in networks
For the case of a network without full permutation symmetry, we consider the following equations of motion for the oscillators θi = σ N j=1 [A] i,j sin(θ j − θ i − α), i = 1, . . ., N, where σ is the overall coupling strength.The adjacency matrix A represents a network that does not have any non-trivial symmetry.To generate this network we initiate the globally coupled network, as depicted in Figure 8 (a), and remove uniformly at random one edge from the graph.The resulting graph has two new clusters that are depicted in different colors, which is calculated via Sage [41].See Figure 9 for an illustration of the (random) network used in our numerical results.We consider time interval of 5000 time steps and sample the trajectory at every 0.02 time steps.Also, the parameters are σ = 0.01 and as before α = π/2 − β with β = 0.025.Choosing the same two-dimensional section as previously, we construct the basin structure via VPS, as illustrated in Figure 9 (b), using five clusters that are determined by the elbow method discussed earlier.Here, we consider a transient time of 4000 time steps and calculate the similarity measure and cross correlation, over the remaining 1000 time steps.
The basin structure shows that the system has a high sensitivity to initial conditions, so our method is able to characterize such basin which has an intermingled structure.It remains an interesting line of research to describe the mechanism behind the chaotic dynamics, since it differs from the globally coupled scenario where a low-dimensional description is available [8].

VPS for phase oscillators
Different from the HR system, when we consider coupled oscillators, the time series is given by phases {θ i (t)} t≥0 in the circle.To construct the VPS, we perform an embedding into R, so we can use the same similarity measure introduced in Section 4.1.More precisely, for each oscillator time series, we consider x i (t) = cos(θ i (t)), i = 1, . . ., N.
To allow the system to settle on to the attractor we use 20, 000 time steps.The final 1, 000 time steps are utilized for the construction of the VPS.Following several of the above examples, we use only the x-component for the VPS.The basin shown contains 1, 000, 000 initial conditions sampled in a grid of [−3, 3] 2 for the x and y component of a randomly chosen node.The initial condition for all other nodes is set to (0, 0) T .Though we are using a different network structure, we find a very similar basin structure to that of [36].This basin appears to be riddled, which indicates how general this phenomenon is.Additionally it is clear that non-trivial symmetry is not necessary for riddled basins.
Figure 10: Basin structure in the Hénon system for a randomly chosen x,y plane.The system is given by Eq. 24.The network is the same DTI network shown in Fig. 1.This structure appears to have a riddled basin.

Figure 2 :
Figure 2: All VPSs in Fig.3 of the main text.(a)-(h) show all e l vectors (Eq.13) in each of 8 clusters.The x-axis represents vertex pairs while the y-axis contains all initial conditions forming forming a cluster.The color intensity represents entries of e l .(i) The color-coded network shows the chimera state in (f).The vertices with the highest degree of synchronization are colored white while the remaining ones are shown in black.(j) shows x i versus time for this chimera state and is plotted as follows: first, x i (t) and x j (t − τ * i,j ) for nodes i and j having the smallest L(i, j, τ * i,j ) value is plotted; thereafter, the phase-shift for any new node k is decided such that L(i, k, τ * i,k ) minimum over all the pre-existing nodes i.

Figure 3 :
Figure 3: Time series of the x-coordinate of the HR neuron model evolving on the DTI network without phase shift.(a)An x-coupled HR neuron model with no chemical coupling (i.e.α = 0 and σ = 1.7),where we only observe identically synchronized oscillators.(b) Zoom in panel of (a) for better visualization.(c) For chemical coupling α = 0.03 (and σ = 1.7) oscillators still spike at the same time, but have different sized peaks creating a spread in the heights, though they still fully synchronize in phase.(d) Zoom in panel of (c) for better visualization.(e)The HR model with the same coupling parameters as (c) but with the model parameter a drawn from a normal distribution a ∼ N (1, 0.1).(f) Zoom in panel of (e) for better visualization.In this case, the oscillators peak in the same general time frame but these spikes have different peaks and different shapes.This represents a more generalized synchronization pattern.

Figure 4 :
Figure 4: (a) The boundary set obtained after simulation of the HR dynamics under coupling given by Eq. 20 which corresponds to the boundary set from Fig.3(b) in the main text.(b) The graph of ln( 1 ϵ ) vs ln N (ϵ), calculating the dimension of the boundary set in (a).The fitted slope of 1.27 estimates the fractal dimension.

Figure 6 :
Figure 6: The elbow method.(a) The traditional version of the elbow method relies on examining the plot of W ( Ĉ) and finding a point discontinuity, which is then chosen as k.Here it is challenging to find a clear elbow.(b) Instead, ln(k) vs. ln(W ( Ĉ)) is plotted, with a linear fit determined.The best k value is chosen to be the last value before which the majority of the points fall on the line of best fit, in this example k = 5.

Figure 8 :
Figure 8: Basin structure of two globally coupled populations.(a) Two globally coupled populations of five nodes.The (intra-)coupling strength inside a population and (inter-)coupling strength between populations are given by µ and ν, respectively, with µ > ν.(b) Two-dimensional section of the state space showing basins of the 18 (clustered) distinct states identified by the VPS of the two globally coupled populations.There is a symmetry among the basins with respect to reflections across the diagonals, which originates from the reflection symmetry of the network.To construct the VPS, we use β = 1 in Equation (13) and a grid of 1248 × 1248 uniformly sampled initial conditions.

Figure 9 :
Figure 9: (a) Network with one single link removed, breaking the full permutation symmetry of the globally coupled network.Nodes with different colors correspond to the new clusters formed from the permutation symmetry group.Red edge corresponds to a randomly selected edge that was removed from the globally coupled network.(b) Two-dimensional section of the state space showing basins of 12 (clustered) distinct states.To construct the VPS, we use β = 1 in Equation (13) and a grid with 1248×1248 uniformly sampled initial conditions.