{'fraction_non_alphanumeric': 0.030386407455012854, 'fraction_numerical': 0.008696176092544987, 'mean_word_length': 4.771763069433175, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Synchronization and clustering in complex quadratic networks

In continuation of prior work, we investigate ties between a network's connectivity and ensemble dynamics. This relationship is notoriously difficult to approach mathematically in natural, complex networks. In our work, we aim to understand it in a canonical framework, using complex quadratic node dynamics, coupled in networks which we call complex quadratic networks (CQNs). After previously defining extensions of the Mandelbrot and Julia sets for networks, we currently focus on the behavior of the node-wise projections of these sets, and on defining and analyzing the phenomena of node clustering and synchronization. We investigate the mechanisms that lead to nodes exhibiting identical or different Mandelbrot set. We propose that clustering is strongly determined by the network connectivity patterns, with the geometry of these clusters further controlled by the connection weights. We then illustrate the concept of synchronization in an existing set of whole brain, tractography-based networks obtained from 197 human subjects using diffusion tensor imaging. Synchronization and clustering are well-studied in the context of networks of oscillators, such as neural networks. Understanding the similarities to how these concepts apply to CQNs contributes to our understanding of universal principles in dynamic networks, and may help extend theoretical results to natural, complex systems.


Background on dynamic networks
Understanding dynamics in coupled networks has become recently an important research direction in a variety of fields. Mathematically, it is a rich topic, that generates interesting questions at the intersection of graph theory and dynamical systems, the complexity of which depends equally on the type of node-wise dynamics, and on the network size and architecture. Mathematical inquiries and results on dynamic networks are relevant to many questions in the life sciences, since many natural systems are organized as complex, adaptable networks. Understanding how the hard-wiring of these systems affects their function is a crucial, structural question, whether the system in question is the brain, a social or an epidemic network.
An important direction of study in dynamic networks relates to their long-term ensemble behavior, and the dependence of this behavior on the system's parameters. These parameters may characterize the individual dynamics of the network's nodes or the connectivity between them, including both connection distribution and connection strengths. For example, if the system describes firing rates in the neural populations that form a brain circuit, one aspect crucial to understand is whether the system converges in the long term to sustained firing (stable equilibrium), or to firing rate oscillations (stable limit cycle). This is important from an application standpoint, since these two different asymptotic patterns reflect into different functional and cognitive outcomes for the circuit in question. Another aspect which was proved to be crucial to a circuit's functional dynamics is that of synchronization in the nodes' activity. Neural populations have been observed to act completely synchronously or asynchronously, or to separate into clusters with different activity patters. On one hand, synchronization and clustering are highly dependent on the network interconnectivity and architecture; on the other hand, at a functional level, they have been shown to reflect onto the type and efficiency of network function to be performed.
These different possible outcomes in dynamic activity of a network depend on intrinsic regulatory parameters that govern the behavior of each unit, in absence of the coupling (node-wise dynamics); but they also depend on the coupling itself: which nodes are connected (e.g. via synaptic pathways in a neural network), and how strong these pathways are in time. Studying this relationship between the network structure and the emergent coupled dynamics is central to the study of dynamic networks. For example, in neuroscience, some psychiatric disorders have been correlated both with atypical brain connectivity patters between specific brain areas, as well as with abnormal neural rhythms in certain these areas. It is important to understand how the former contribute to the latter, and how this combination further leads to a certain ensemble behavior in the brain network versus another.
This tie between a network's hard-wiring and function has been studied intensely for the past decade (e.g., in neuroscience, through the Human Connectome Project). It has become increasingly clear that establishing a precise correspondence between connectivity and emergent dynamic patterns is a very difficult problem both analytically and computationally, especially when considering the sheer size and complexity of natural networks, with the brain in particular. Newer research has instead focused on understanding the principles of this correspondence in smaller, simpler networks (in terms of network size and architecture, or in terms of node-wise dynamics, or both).

Complex quadratic networks, as an analytic framework
In our previous work, we adopted the simplified framework of networked logistic maps as a starting point for approaching basic dynamic questions in the context of networks [1,2]. That is, we considered nodes with complex quadratic dynamics in the family f c : C → C, f c (z) = z 2 + c. More specifically, in this framework, each network node receives weighted inputs from the adjacent nodes, and integrates these inputs in discrete time as a complex quadratic map. Then the system takes the form of an iteration in C n : where n is the size of the network, A = (A jk ) n j,k=1 is the binary adjacency matrix of the oriented underlying graph, that is A jk = 1 if there is an edge from the node k to the node j, and A jk = 0 otherwise. The coefficients g jk are the signed weights along the adjacency edges (in particular, g jk = 0 if there is no edge connecting k to j, that is if A jk = 0). In isolation, each node z j (t) → z j (t + 1), 1 ≤ j ≤ n, iterates as a quadratic function f j (z) = z 2 + c j . When coupled as a network with adjacency A, each node will act as a quadratic modulation on the sum of the inputs received along the incoming edges (as specified by the values of A jk , for 1 ≤ k ≤ n).
This is a simplified framework, which turned out to preserve, often in a weaker form, some of the properties and results determined in the traditional case of a single iterated quadratic map (e.g., the existence of an escape radius is guaranteed in some types of networks, but may fail in others). We used this to our advantage, and proceeded to study theoretically the effects of network architecture on its long-term dynamics. We have been generally investigating whether one can use properties of multi-dimensional orbits in C n , in particular their asymptotic behavior (via the topological and fractal structure of Julia and Mandelbrot multi-sets), to classify dynamic behavior for different network architectures. To this aim, we defined network asymptotic sets as follows: Definition 1.1. We define the multi-Mandelbrot set (or the multi-M set) of the network the parameter locus of (c 1 , ..., c n ) ∈ C n for which the multi-orbit of the critical point (0, ..., 0) is bounded in C n . We call the equi-Mandelbrot set (or the equi-M set) of the network, the locus of c ∈ C for which the critical multi-orbit is bounded for equi-parameter (c 1 , c 2 , ...c n ) = (c, c, ...c) ∈ C n .
In our previous work, we focused on describing topological properties of these asymptotic, global network sets. We aimed to understand which geometric aspects survive when generalizing these objects from the context of a single iterated map to the richer dynamics of ensemble iterations of a whole network.The difficulty comes from the fact that, even in the situation where the nodes are identical, the connectome favors nodes differently -hence the strong dependence of the topology of the asymptotic sets on the architecture of the network. We have been therefore trying to track down how the connectivity profile of the network affects the shapes, properties and relationships between asymptotic sets.
For example, note that the Mandelbrot set of an n-dimensional network was defined as the node parameter range for which the critical point (i.e., all nodes equal zero) is bounded in C n under ensemble iterations of the network. In the traditional case of single iterated quadratic maps, this is equivalent to Julia set connectedness locus. One of our central questions in previous work has been whether this result remains true for networks. Numerical experiments strongly suggest that this equivalence no longer holds in the case of iterating networks (for either Julia sets, or uni-Julia sets). In our previous work, we conjectured a weaker form of this result, for the case of identical nodes. We suspect (but have not yet proven analytically) that connected uni-J sets only occur for c in the network equi-M set, and that totally disconnected uni-J imply a c outside of the equi-M set. These seems to be a transitional region c ∈ C for which the uni-J set is disconnected, but not totally disconncted; this region contains the border of the network's equi-M set. This is not surprising, since in order for a network's critical orbit to be bounded, all of its node-wise components need to be simultaneously bounded, a complex condition that depends crucially on the shape of the network. A natural question that emerges is that of how different can node-wise behaviors be, between different nodes of the same network. We have shown before that even in networks with identical node dynamics (i.e., all nodes have the same parameter c) the nodes may exhibit different asymptotic behaviors (in the sense of their own, node-wise, Mandelbrot and Julia sets). Node coupling determines the communication between two or more nodes, and different types of networking motifs seem to lead to network Mandelbrot and Julia sets of different shapes and sizes. In the current paper, our interest is focused on understanding node "synchronization." We aim to analyze and interpret the distinct effects of varying connection patterns on synchronization. We are also interested to understand the relationship between node de-synchronization, and failure of the Fatou-Julia theorem.
2 Node-wise synchronization and clusters in CQNs

Definitions
In all the following sections we will be working with n-dimensional networks with complex nodes z 1 · · · z n with identical node-wise dynamics (i.e., all nodes are governed by the same complex parameter c).
Definition 2.1. For each node z k , 1 ≤ k ≤ n, we define the node-wise equi-M set M k as the parameter locus c ∈ C for which the component corresponding to the node z k of the critical multi-orbit is bounded.
Note that, with this definition, the network equi-M set is the intersection of all node-wise equi-M sets. Part of our mission in this paper is to understand what conditions may determine whether the equi-M sets for two nodes are identical or different. To help with this aim, we introduce some additional terminology: Definition 2.2. We say that two network nodes z i and z j are M-synchronized, if their equi-M sets M i and M j are identical. We also say that the nodes belong to the same M-synchronization cluster (or simply M-cluster) of the network.
and all unit weights g ij on the nonzero edges (i.e., the adjacency matrix is also the connectivity matrix). In other words, the corresponding 5-dimensional discrete iteration is given by following scheme (where all nodes update simultaneously): The five nodes separate into three different M-clusters, the M sets of which are shown as contours in Figure 1b: z 1 , which is only driven by self-input, generates a node-wise M set identical with the traditional Mandelbrot set (the boundary of which is shown in red). Nodes z 2 , z 4 and z 5 are M-synchronized into one M-cluster; their identical M sets are sketched with a blue contour. Node z 3 forms the third M-cluster, the set of which is shown as a green contour.
The concept of synchronization has been widely studied in networks of continuous-time dynamic oscillators, with brain networks in particular. In this more traditional form, synchronized activity of two nodes is understood as convergence to an identical temporal oscillatory pattern. To better understand how our concept of M-synchronization relates with the standard definition, we illustrate the latter in the context of the Kuramoto model (as an example of a simple, generic model of network dynamics). Synchronization patterns in Kuramoto coupled dynamics, and how they are actualized in the brain connectome have been studied over the past decade [3]). For similarity with CQNs with identical nodes, we consider the Kuramoto model with N identical coupled oscillators where the variables θ i can be viewed as angles modulo 2π. Here, w is the common internal frequency, K is the control parameter, and A represents (similarly to our model) the connectivity weight matrix. Depending on the context of w, K and A, subsets of nodes in the network may converge in time to identical or out of phase oscillations in [0, 2π], creating complex synchronization patterns. The N -dimensional network can be then subdivided into subsets (or clusters) of temporally synchronized nodes, ranging from one cluster (in the case where all nodes are eventually synchronized) to N distinct clusters (when all nodes perform, in the long term, distinct/out of phase oscillations). As an initial comparison, we show in Figure 1c an example of clustering in the Kuramoto network with the same connectivity matrix as the one used for the M-clusters.
In this particular example, the activity in the Kuramoto network readily synchronizes into the same three clusters as the CQN, as shown. However, this cannot be expected to happen in general; in fact for smaller values of the coupling K, the system may undergo other transients along the way to this asymptotic pattern, making it hard to even distinguish when the oscillations have settled around the stable periodic attractor. We will return to this phenomenon and its significance later in the paper.

Effects of networking on M-clusters
We have previously shown that topological properties of the network equi-M set in CQNs depend strongly on the network connectivity patterns. For example, we have illustrated that, even in low-dimensional networks, small perturbations in connectivity (understood as small local changes in either position of the connecting edges, or in the connectivity weights) may lead to large differences in the size, shape and connectedness of the equi-M set. Here, we wish to establish that these differences stem in fact from the clustering behavior. To begin with, it is easy to see from our definitions that the equi-M set of a network, is the (finite) intersection of all its node-wise Mandelbrot sets. This relationship may be even more meaningful than seems at first. For example, it appears from our simulations that, for many networks, this is a nested intersection. Below, we describe a 3-dimensional example, which we will also revisit later in the paper. Example 1. In the network defined in Figure 2, the nodes define three distinct, nested equi-M sets.
Proof. We show first that, if the critical orbit is bounded in the z 1 component, it is also bounded in the z 3 component. Indeed, for any iteration n ∈ N, we have the following: Since z 1 is bounded, it follows that z 3 is also bounded. Consider now the critical orbit of the network under the parameter c = −2. It is easy to show that this critical orbit is eventually fixed in its z 3 component (z 3 (n) = 2, for all n ≥ 2), hence −2 ∈ M(z 3 ). We will show that the z 1 component of the critical orbit escapes. Notice first that the sequence z 1 (n) has all positive real values, and that, for n ≥ 2 Since z 1 (n + 1) z 1 (n) > 2 for n ≥ 2, it follows that z 1 (n) → ∞. This proves that M(z 1 ) M(z 3 ). For the other connected node pair, one can show that M(z 2 ) M(z 1 ) very similarly to the inclusion above. In addition, c = −3/4 ∈ M(z 2 ), but 3/4 ∈ M(z 1 ). To show this, notice that the z 3 component of the critical orbit is eventually fixed to −3/4, and the z 2 component of the critical orbit eventually oscillates between −3/4 and 3/2. Hence the z 1 component is a sequence of positive real numbers, which builds up recursively, for every even n ∈ N, as: which clearly diverges.
The nesting property does not hold in general. One can simply think of the trivial example in which the nodes are independent, with different levels of self-input: Notice that the components of the critical orbit scale linearly with the value of the parameter; hence the node-wise M sets also scale: if c ∈ M(z 1 ), then a1 a2 c ∈ M(z 2 ). Then, for values of a 1 and a 2 close enough together, the sets intersect nontrivially.
There are, however, more complex examples of networks for which there exist two nodes with M sets intersecting non-trivially, as shown below. Figure 3 illustrates the geometry of the M set clusters for the N = 20 dimensional network with adjacency matrix A1 (shown in Appendix B), and unit edge weights. The network exhibits eight distinct M clusters (shown on the left), at least two of which (illustrated separately on the right, for clarity) intersect non-trivially. In future work, we are interested to further investigate necessary and sufficient network conditions for nesting of node-wise uni-J and M sets. In the 20-node network with adjacency matrix A1 (shown explicitly in Appendix B) and unit edge weights, the nodes synchronize into eight distinct M clusters (containing from one to four nodes each). Their contours are illustrated in different colors. The contours are not all nested, i.e. some of them intersect non-trivially. Right. The panel shows in more detail two of the M cluster contours on the left which have a non-trivial intersection; we specifically marked one point (c=−0.23 + 0.05i, in red) which is in the blue cluster, but not in the maroon cluster, and one point (c=−0.2 + 0.06i, in green) which is in the maroon cluster, but not in the blue one. The computations supporting this panel can be found in Appendix A.
Additional robust properties of M sets have been previously pointed out, such as the persistence of the cusp, despite the fact that network M sets don't otherwise preserve the hyperbolic bulb structure of the traditional Mandelbrot set. While the shape and structure of the tail varies wildly with that of the network (including the potential to break off into separate connected components), the cusp is a common geometric feature in every node-wise M set.
In this context, it becomes clear that, to understand how the geometry of the network equi-M set depends on connectivity patterns, it is helpful to first understand how the number and shape of the node-wise clusters depend on these patterns. In particular, we are interested to study the differences between the effects on clustering of perturbing global connectivity (number of edges), versus local connectivity patterns (removing, adding or replacing one edge), and versus altering the connection strengths (increasing and decreasing one or multiple edge weights).
Local changes in edge placement and weights. Our simulations suggest that simply strengthening or weakening one edge can alter the shape of the node-wise M-sets in one or more clusters, while leaving the cluster separation unchanged. In contrast, adding or removing one edge completely can promptly change the cluster structure, in addition to potential changes in the shape of the sets corresponding to each cluster.
To illustrate, consider the 5-dimensional matrix obtained from the connectivity matrix A described in example (1), to which we add cross-talk from the node z 3 to the node z 1 . The panels of Figure 4 show how the shape of the clusters changes when the weight of this new edge decreases from A 13 = 1 to A 13 → 0. Notice that, throughout this process, the two-cluster structure remains unaltered, although the shape of one of the clusters is affected by the perturbations in the local weight. As A 13 becomes zero, a bifurcation occurs in the clustering structure, to recover the original 3-cluster formula in Figure ?? (due to the fact that node z 1 has now been decoupled from any input coming from other network nodes and is only being driven by the self-loop, hence the associated M set is the traditional Mandelbrot set). In Figure 5, we continue by illustrating how local connectivity changes (consisting of removing or adding one edge), applied to the matrix A, affect clustering and the shape of the M set in each node. For example, removing the self-loop on node z 3 (i.e., setting A 33 = 0) allows the node to synchronize with the others in the second cluster for the original matrix A, and the resulting network will now only have two clusters (panel A). The shape of the M set for this enlarged cluster is different, however, than the shapes of either the two merging clusters for the original matrix A. Adding an equal weight cross-talk from node z 3 to z 2 (i.e. setting A 21 = 1) has the same effect of synchronizing nodes z 2 − z 5 , into a cluster with a new shape for the M set. Finally, replacing the self-input on node z 1 by a same weight input from node z 2 synchronizes z 1 with the rest of the nodes in cluster 2 of A, without changing the shape of the remaining clusters.
Continuing the comparison in Section 2.1, Figure 6 shows the synchronization clusters formed in the Kuramoto model with five identical nodes and the same connectivity matrices as in Figure 8. The clustering patterns in the Kuramoto model depend on the initial conditions and on how these fall into the attraction basins of various stable oscillations. In particular, identical initial conditions produce trivial clustering behavior (nodes immediately synchronize to one cluster, independently on the network). By analogy with the criticality condition in equi-M sets, we considered in for our simulation in Figure 6 initial conditions with balanced phases around the circle, for which the order parameter is zero. Even with fixed initial states, the temporal evolution of the Kuramoto system is very complex, and quite difficult to interpret. Even for the simple networks at hand, the system exhibits chimeras, with synchronization patterns evolving in time, making it difficult to pin down asymptotic states, and even more difficult to compare the system's long-term behavior between different connectivity schemes. Figure 6 illustrates this comparison for the same three matrices used in Figure 8 for CQNs. Notice that, since they start with the same initial conditions, the oscillations are rather similar between the three examples for the first transient segment (panels 1). However, Figure 5: Effect of edge position on equi-M set clustering and shape. A. Deleting the input A32 = 1 leads to two clusters: cluster 1 (node 1, in blue) and cluster 2 (all other nodes, in red). B. Adding A33 = 1 leads to the same two clusters. C. Deleting the input A11 = 1 and adding the input A12 = 1 leads to two different clusters: cluster 1 (nodes 1,2,4,5, shown in blue) and cluster 2 (node 3, shown in red).
these oscillations weave into different evolving synchronization patterns, according to the underlying matrix. For the first connectivity scheme A, nodes 3,4 and 5 initially synchronize (panel 2), then de-synchronize (panel 2), then nodes 1,2 and 4 synchronize, showing three clusters at the end of our simulation. Even though the final result is similar for the second connectivity scheme (three clusters and the same grouping as above), this system crosses a 4-cluster transient. Comparing the two systems prematurely in time may suggest different behaviors, when at the end of the simulation they are very similar. Conversely, our third example shows four transient clusters around t=250, same as example 2, but its asymptotic behavior is completely different (all nodes fully synchronize by the end of the simulation run t = 1450. In conclusion, while it is clear that the connectivity matrix A controls to a large extent the Kuramoto patterns, the assessment of how and to what extent this is the case relies on choosing the appropriate moment in time to make this assessment, which is an additional complication to an already complex problem. When working with CQN networks, a good approximation of asymptotic behavior is guaranteed in the case of equi-Mandelbrot sets by escape radius theorems, eliminating in most cases the question of whether the set would further change under prolonged iterations. How the behavior of the balanced initial conditions is relevant to synchronization in the Kuramoto model for more general solutions is a well-studied area, beyond the scope of our paper. However, the analogous question in CQNs, asking to what extent the critical orbit truly encapsulates the behavior of all possible network orbits, is a new and important extension of traditional results in single quadratic map iterations, and is further discussed in Sections ??. A more careful comparison of parameter and network dependence between the two models (CQN and Kuramoto) is carried out in Section 2.4. , for 1 ≤ i ≤ 5, as the connectivity matrix undergoes the same changes as those described in Figure 5. In order to better emphasize how synchronization patterns change over time, we selected to illustrate specific time windows, separated by grey shaded blocks. Global changes to edge density. On a more global scale, our simulations suggest that the number of edges correlates with both the size of the M sets and their clustering, in the sense that a higher number of edges leads to smaller and more synchronized node-wise M sets. This is not surprising, and has been suggested in our prior work. We illustrate this relationship further in Figure 7, in which we consider the behavior in three different networks with N = 20 nodes, unit edge weights, and progressively higher edge densities (from panel a to panel c). The network in panel a shows eight different clusters, with at least two M sets intersecting non-trivially. The network in panel b exhibits three nested M sets, and panel c illustrates a network where all M sets are synchronized. This follows the intuition that, the denser the matrix -the more synchronized the M sets. We will reiterate this idea when illustrating synchronization with an example from brain connectomics, in Section ??.
Of course, a more comprehensive analysis needs to be performed to establish whether these effects remain significant across all possible matrices with these properties. For the start, this is more easily done for lower dimensional network spaces, where there are less configurations to investigate. In the next section, we consider all 3-dimensional oriented networks with self loops, and fixed edge weights (=1), and we study more carefully the effect that the re-positioning of any subset of edges produces on the cluster structure and on the geometry of the corresponding M sets. Figure 8 illustrates the shapes of the node-wise M sets and their clustering behavior for all possible oriented networks of three nodes, all three of which have self-loops, and such that all active edge weights = 1. Via symmetries, there are 16 such networks. Each figure panel corresponds to one of these networks (shown in the top left corner without the self loops, to avoid cluttering). We will refer to these panels in sequential order 1-16, from left to right and then top to bottom.

Classifying synchronization in 3-dimensional CQNs
For each network, the system's M-set may synchronize into one, two or three M-clusters. We used different colors (blue, red and green) to distinguish the M-sets in different clusters, as well as to label the network nodes that correspond to each individual cluster. The first node z 1 is always labeled as "cluster 1" and marked in blue, together with all nodes that synchronize with it. The first node (in increasing order of the indices) that is not part of cluster 1 defines "cluster 2" and is marked in red, together with all the nodes that synchronize with it. If there is any node left that is not already in clusters 1 or 2, it automatically forms cluster 3, and is marked in green. Seven of the 16 cases lead to a single M-cluster, seven of them lead to two M-clusters, and the remaining two lead to the formation of three M-clusters. In some cases, the mechanisms behind the clustering behavior are intuitively clear; in other cases, the explanations are more subtle, or even counter-intuitive. For example notice, rather trivially, that symmetries between nodes lead to synchronization of those nodes into the same cluster. For example, in panels 1, 3 and 6, all nodes are identically coupled in the network, hence they have identical dynamics and they all synchronize into one cluster. It was pointed out in prior studies that higher node degrees lead to smaller network M sets, which is also reflected here. In panels 2, 9, 10, 13, it is also easy to notice why two of the nodes have identical wiring, hence they synchronize, while the behavior of the third node is different, and defines a second cluster. It is still somewhat intuitive why there are three different clusters for a network like that in panel 8 (representing a feed-forward network with self loops): node z 3 is independent on the other nodes dynamics, hence its M set is the traditional Mandelbrot set; z 1 feeds directly on z 3 , and the integrated self and external input lead to a smaller M set for node z 1 ; in turn, z 1 projects onto z 2 , leading to an M set for z 2 that is a subset of that corresponding to z 1 (hence the three clusters).
In this simplified scheme that analyzes only three-node interactions, it is therefore more tractable to tease apart which aspects of the network wiring lead to which which clustering scheme. Once the cluster formation is justified in each case, one interesting point of discussion is to track how changes in networking trigger synchronization changes. For example, if the connection from z 1 to z 2 in panel 8 is used instead to connect z 3 to z 2 (panel 9), then the node z 2 will no longer form its own cluster (with smaller M set than that of z 1 ), but instead it will join z 1 in its cluster, with no modifications to its M set. Interestingly, the same exact effect is produced on the M sets if the edge is simply removed altogether (panel 7). If, on the other hand, one restarts with the network in panel 8 and adds feedback responses along both inter-node inputs, all nodes synchronize to a new, smaller M set (panel 5). If one adds feedback just along one of these two edges, the network will switch to a two cluster structure (panel 12).

Comparison with the Kuramoto model
In Section 2.2, we illustrated a simple example of how perturbations in connectivity affect synchronization patterns in a Kuramoto network. Here, we focus on describing similarities and differences between synchronization in Kuramoto networks and equi-M set synchronization in CQN networks. We discuss possible advantages of using equi-M set synchronization as a classification tool over the standard synchronization techniques.
Instead of embracing standard synchronization in temporal activity, in CQNs we use equi M sets to describe synchronization of the nodes' critical behavior (bounded versus unbounded) in the complex parameter plane. For simple iterations of polynomials, it is well known that the critical orbits encode information about all accessible orbit behaviors. For example, for traditional iterations of quadratic functions in the family f c (z) = z 2 + c, the Fatou-Julia theorem states that, for values of c in the Mandelbrot set (i.e., the iteration of the critical point ξ 0 = 0 is bounded) all other initial conditions with bounded orbits ("prisoners") form a connected set around ξ 0 . In previous work, we have started exploring to what extent such a characterization remains true in the case of networks [2]. In the meantime, equi-M sets represent a way to visualize the sets of parameters for which different CQN nodes have different behavior.
Just like for CQNs with identical nodes, where the node dynamics are specified by the complex parameter c = (Re(c), Im(c)), the Kuramoto network with identical nodes (2) is controlled via a two-dimensional parameter plane (w, K). By analogy with the classification of M-clusters in 3-dimensional CQNs, shown in Figure 8, in Figure 9 we illustrate the parameter plane classification of phase coherence for the corresponding 3-dimensional networks of Kuramoto oscillators. More precisely, for each of the same 16 adjacency matrices described in the case of CQNs (Figure 8), for each parameter pair in the region (w, K) ∈ [0.1, 1.1] × [0.1, 1.1], we considered the orbit initiated at balanced initial conditions (2π/3, 4π/3, 2π), and we estimated the number of synchronization clusters after T = 500 seconds. Notice that, as previously discussed, these results may change when T is increased, in ways which depend on the parameter values and on the network itself (as further illustrated in Figure 10). This all suggests that our M-synchronization approach is a more efficient tool when evaluating how emerging dynamic behavior depends on the network structure. For example, there are eight identical "phase coherence" parameter planes in the Kuramoto model simulations, all with R = 1 (identifying a unique synchronization cluster throughout the region of interest in the (w, K) plane). These correspond to widely different connectivity matrices, with different graph theoretical properties, differences which may be relevant to brain function and behavior, and which are not picked up by the synchronization patterns in a network model with Kuramoto node dynamics. As a very particular example consider networks 8 and 11 in the grid in Figure 8 (read top-down, then left to right). They represent a two-link motif and a triangular motif, respectively; the difference between the strengths of these motifs are believed to be a potentially important factor contributing to modulating brain activity and behavior. The Kuramoto model does not distinguish well between the patterns generated in these two networks. In turn, the M-synchronization technique in CQNs found an interesting difference. Even though the system has three M-clusters in both situations, the equi-M set corresponding to the node with the additional input in the triangular motif is small than in the two-link motif (a subset). In empirical work, we identified the size of the equi-M set to be a marker for gender and other physiological differences. It is very promising that such a marker can be tracked down to the shape of the network when using CQNs. Among other factors, one advantage of CQNs consists of being able to observe the behavior of each identical node (critically bounded or not bounded) independently from the others, without ignoring the connectivity web between them. Moreover, the geometric properties of equi-M sets represent as many ways to assess the effects of the underlying network on dynamics -a set of measures which is substantially more complex and subtle than the phase coherence tools in Kuramoto networks. While the two models will always describe different types of dynamics, the M-synchronization appears to be a model that is better equipped to answer questions on the relationship between network architecture and the emerging coupled dynamics.  For all simulations, the intrinsic frequency w was fixed to w = 1.

M-synchronization in tractography-based networks
To bridge this concept with its potential practical applications to neuroscience, we use human tractography data to illustrate a practical example where M set synchronization can be used to understand the role of weak connections within the brain connectome.
Data description and processing. Our analysis is based on tractography-derived, neural connectivity data obtained from the S1200 release of the Human Connectome Project [4]. The data set and the structural network construction are described in detail in Appendix C. In resulting structural networks provide us for the mathematical analysis with non-oriented, unsigned, normalized connectomes for 197 subjects. For each subject, the brain was partitioned into 116 areas (which we use a our network nodes), and normalized connection weights were obtained between each node pair (as described in Appendix C). Since these connectomes are tractography-based only, they do not capture either the directional aspect or the signature (excitatory vs. inhibitory) of the information transfer between node pairs. The corresponding matrices are therefore symmetric matrices with positive entries. For this mathematical exploration, we chose to "enhance" the connectome of each subject by adding the identity matrix to the tractography-based connectivity matrix (i.e., we filled the first diagonal with unit entries). We did so for two reasons. First (a mathematical reason): the inputs coming from other nodes are collectively very weak, preventing many receiving nodes from escaping to infinity, and thus leading to many trivial node-wise M sets (equivalent to the whole complex plane). Second (a modeling reason): this also reflects the corresponding situation in real brain networks, in which nodes self-regulate, even though these "self-loops" cannot be captured by tractography data.
From each of these "enhanced" matrices, we used Matlab to compute the node-wise M sets for all the 116 nodes (using 50 iterations of the network, with spacial resolution 200 × 200 pixels in the c-plane). We identified two nodes as being M-synchronized based on their M-sets occupying the same locus (via a maximum error = 10 pixels). Using the original connections, all the nodes were synchronized (i.e., showing only one M set cluster). However, when progressively pruning out the weakest connections, we observed the nodes subsequently de-synchronising (see Figure 11 for an illustration). This is not surprising, since a low level of cross-talk between nodes is expected to contribute to their synchronization, hence eliminating some of that cross-talk would lead to a more diverse synchronization profile. However, as one would suspect, this increasing relationship between the number of clusters versus the threshold value cannot possibly extend universally, since eventually some nodes will be stripped of all interconnections, and, being driven solely by the unit self-input, will synchronize to the traditional Mandelbrot set.
To more consistently follow the effect of eliminating weak connections on synchronization, we used a moving threshold h which we gradually increased from h = 0.01 to h = 0.1 in increments of ∆h = 0.005. For each value of h, we then counted the number of different node-wise equi-M clusters, when deleting all connections weaker than the set threshold h. To give an idea of the effect of such thresholding on the connectome, let us notice that the low level cross-talk is prevalent. While the maximum strength connection among individuals ranges from 0.15 − 1.04, the mean weight of active connections among individuals ranges between 0.006 − 0.012. Hence when we delete, for example, all connections weaker than h = 0.01 , only the strongest 16-32% connections are retained (with different individuals placing differently within this range). When deleting all connections weaker than h = 0.03, we are left with only the top 4-11% connections; and when deleting all connections weaker than h = 0.05, only the strongest 1-5% connection are left.  With this in mind, we see the progression of de-synchronization as the threshold h is being increased. Figure 12 shows the mean number of distinct M clusters across the subject population, with error bars showing between-subject variability, and how this mean evolves with respect to the increasing threshold h. Notice that the maximum desynchronization is reflected as 11-23 clusters (with variability between individuals, as shown in the histogram in Figure 13a), and these maxima are achieved for more than half of the subjects within the interval h = 0.05 − 0.06 (the variability between subjects is shown in the histogram in Figure 13b).
First, the maximum desyncrhonization is consistently achieved when most of the cross-talk between nodes has been eliminated, and only the strongest 1-5% of connections are left. While this may be initially counterintuitive, it in fact relates to existing knowledge of information processing in dynamic networks. The idea that sparse networks are useful for brain processing, and may in fact outperform denser networks in certain contexts is not new [5]. In sparse networks, the contributions of individual nodes are less centralized, allowing more flexibility of the ensemble response. In our context, this reflected in the fact that sparser networks exhibit a larger number of clusters, while denser networks are more likely to lead to node synchronization, and diminishing the within-network variability that may be required for accomplishing certain brain functions.
The low-level cross-talk between brain areas may allow the brain to swiftly switch between synchronization modes, in order to produce a desired functional outcome. Our results suggest that the physical mechanism used to sustain this possibility may involve maintaining a multitude of weak connections, which can be used as needed for active information transfer, based on context and desired outcome. Broad and quick decisions on how much cross-talk to actively involve in a certain process is consistent with our knowledge of both brain plasticity and dynamics. This mechanism may be a cost reasonable choice for the system to efficiently accomplish synchronization switching, since weak connections are more inexpensive to maintain, and easier to sacrifice. To determine which of the anatomical connections between nodes are actually in use at each moment, would require an imaging modality that reflects functional connectivity.

General comments on synchronization and clustering
In this paper, we explored the potential of implementing a synchronization concept driven from applications in network neuroscience, to the theoretical context of networks with complex quadratic nodes. The aim of this is to facilitate analytical approaches to the complex phenomena in nonlinear networks of coupled neural oscillators, by drawing analogies with canonical networks with more tractable dynamics (such as CQNs).
We determined that network nodes separate into a finite number of clusters, based on the criterion that members of each cluster have identical M sets. The cluster number and grouping are primarily determined by the distribution of edges in the network, and the geometry of the M sets in each cluster is further tuned by the connectivity patterns and weights. We found that the network equi-M set (defined as the intersection of the node-wise M sets) may or may not be a nested intersection, depending on the network architecture.
One significant limitation relates to interpreting some of the illustrations based on numerical computations. That is because the existence of escape radius has been determined for a few network types, such as feed forward networks with self loops and dominated networks, but is not known in general. This requires additional caution when interpreting some of the illustrations, which are based on a finite number of iterations; the potential exists for the orbits to remain bounded after escaping the large disc set in the simulation code as the escape domain. Establishing and calculating escape radii for more network classes is one of the primary goals of our future work, and will increase the confidence when interpreting simulations.

Connections with other models and applications
Part of our goal when studying CQNs consists of understanding to what extent methods from traditional, single map discrete dynamical systems can be extended to the context of networks. The other aspect of this approach is that of understanding how behaviors in CQN theoretical networks may be universal, or similar phenomenologically to behaviors of networks of interest in the natural sciences. Translating phenomena such as "node synchronization" and "clustering" into the framework of CQNs is a significant step in this latter direction.
The idea that network connectivity governs the nodes' synchronization and clustering is a central theme in computational neuroscience. For example, modeling work in basal ganglia networks found that different types of connectivity patterns lead to different synchronization clusters and spiking rhythms [6]. Other studies show that even in networks with fixed architecture (all-to-all inhibitory reticulate thalamic networks) small perturbations in synaptic conductance can trigger restructuring of the nodes into different numbers of synchronized clusters [7]. While our model is very different from the ones above, both in the type of dynamics (discrete versus continuous time), and in the definition of synchronization, our results capture efficiently dependencies of clustering on local and global perturbations to the network structure and weights. In fact, a simple comparison suggested that our CQN model outperforms Kuramoto networks in terms of capturing how network perturbations affect synchronization patterns.
Our exploration of the brain tractography data emphasized the idea that weak background cross-talk between various parts of the brain plays an extremely important role in syncrhonizing dynamics. This has been observed before in networks (e.g., opinion formation in social networks), but has never been fully explained. A complex network like the brain has to function within an optimal range, where it can swiftly adapt to "synchronize" or "de-synchronize," and may be able to do so by using small changes to the level of cross-talk. This study shows how our abstract mathematical approach can be used to capture the link between weak connections and brain synchronization patterns.
A future direction, tying into applications and personalized assessments, is to further study betweensubject differences in synchronization profile and threshold effects. As Figures 12 and 13 show, there is variability between the number of synchronization clusters for different subjects for any given threshold. Moreover, the maximum number of clusters differs in value, and is achieved at different thresholds between individuals. This likely reflects different efficiency levels between individuals when switching synchornization regimes, and may be significant to which dynamic patterns are more likely to occur in each subject.
More generally: as we are slowly starting to learn more about network dynamics, reservoir pools, deep learning and artificial intelligence, a lot of important questions in these fields yet remain unanswered by the traditional methods and techniques associated with dynamic networks. Our study illustrates how complex dynamics and our extension of Julia-Mandelbrot theory can be used to contextualize and address some of these questions in a tractable framework.

Appendix C: data and structural network construction
The project released to the public domain extensive MRI, behavioral, and demographic data from a large cohort of individuals (> 1000). For our own study, we considered a subject subgroup of N = 197 individuals (90 males and 107 females, ages from 22 to over 36 years old). The original Structural Network Construction Preprocessed diffusion weighted images were used to construct structural connectomes for the subjects. Fiber tracking was done using DSI Studio with a modified FACT algorithm [8]. As a first step, data were reconstructed using generalized q-sampling imaging (GQI) [9]. Diffusion weighted images were reconstructed in native space and the quantitative anisotropy (QA) for each voxel was computed. Fiber tracking was performed until 250,000 streamlines were reconstructed with angular threshold of 50o, step size of 1.25 mm, minimum length of 10mm, and maximum length of 400mm. Streamline counts were estimated for two parcellations schemes based on the AAL [10] atlas version-1 containing 116. The AAL atlas was registered to the isotropic diffusion component (ISO) image, an output of GQI. Registering directly to the ISO image minimizes any registration issues that could arise by first registering to an individual's T1w image. Atlas registration was conducted using FSL FLIRT [11] with default parameters. A template T1w image was registered to the subject-specific ISO image. The transformation matrix obtained from the registration was applied to all regions in the AAL atlas. In case two regions were registered to the same voxel, the voxel was assigned to the region with the highest probability. In order to fully sample the fiber orientations within a voxel, tracking was repeated 250x with initiation at a random sub-voxel position generating 250 connectivity matrices per subject. For each subject, an undirected weighted structural connectivity matrix, A, was constructed from the connection strength based on the number of streamlines connecting two regions The final connectivity matrix included edges which were present in at least 10% (or 25) of the 250 matrices. Connectivity matrix was normalized by dividing the number of streamlines between each two coupled regions, by the combined volumes of the two regions.