Complex Network Geometry and Frustrated Synchronization

The dynamics of networks of neuronal cultures has been recently shown to be strongly dependent on the network geometry and in particular on their dimensionality. However, this phenomenon has been so far mostly unexplored from the theoretical point of view. Here we reveal the rich interplay between network geometry and synchronization of coupled oscillators in the context of a simplicial complex model of manifolds called Complex Network Manifold. The networks generated by this model combine small world properties (infinite Hausdorff dimension) and a high modular structure with finite and tunable spectral dimension. We show that the networks display frustrated synchronization for a wide range of the coupling strength of the oscillators, and that the synchronization properties are directly affected by the spectral dimension of the network.

Network Theory has provided a large evidence that network topology affects the dynamics and function of complex networks [1][2][3][4][5][6][7] . Recently, growing attention has been devoted to analyze networks not only from the topological viewpoint but, also, from a geometrical perspective [8][9][10][11] . Characterizing the geometry of complex networks is particularly relevant for brain research 12,13 , where the embedding in a three-dimensional space is essential to understand the wiring diagram characterizing the connections between brain regions -or connectome [14][15][16] and also at a more microscopic neuronal level. Interestingly, recent experimental results have shown that the synchronization properties of neuronal cultures grown on 2D slices differ considerably from those grown on 3D scaffolds; these last ones turn out to be much more likely to maintain synchrony and have two regimes: a highly synchronized regime and a moderately synchronized regime 17 . This result can be compared with results obtained in the framework of the Blue Brain project where it has been observed that pairs of neurons have a more significant correlations in their dynamics if they belong to higher dimensional simplicies 18 .
These results reveal that not only network topology but also network geometry plays a crucial role in determining the dynamical properties of a network. However, the connection between network geometry and dynamics has been so far mostly unexplored. Here, we provide a theoretical framework based on simplicial complexes 13,[19][20][21][22][23][24] to investigate the effect of network geometry on synchronization dynamics. Our study reveals that a finite spectral dimension of geometrical networks can be combined with a complex modular network structure to give rise to frustrated synchronization and spatio-temporal fluctuations of the order parameter of the synchronization dynamics. In this way, we provide a clear evidence of the important effect that network geometry has on the dynamics of complex networks.
Synchronization phenomena are the subject of extensive research in physical, biological, chemical and social systems 25,26 . The synchronization properties of a network are influenced by the network topology 27 and by its spectral properties 28,29 . For instance the stability of a fully synchronized state is known to depend crucially on the ratio between the Fiedler eigenvalue (i.e. the smallest non zero eigenvalue) and the largest eigenvalue of Laplacian 28,29 .
Interestingly, complex network topologies characterized by a hierarchical modular structure have been shown to display a dynamical phase called frustrated synchronization 30 in which the order parameter at large times is not stationary but it is instead affected by large spatio-temporal fluctuations. This exotic phase occurs in highly modular networks where contact processes have also a dynamical behavior dominated by rare regions of activity in correspondence with the so-called Griffith phase [31][32][33] . However, in ref. 30 the considered network structures have a finite Hausdorff dimension while a large number of brain networks is known to be small-world 2,12 , (i.e. have an infinite Hausdorff dimension). Therefore exploring whether such dynamical phase can appear in small world networks, such as those that characterize the structure of the brain, has particular relevance in the field of neuroscience and neural computation.
Here, we show that it is possible to generate network geometries called Complex Network Manifolds 19,20,22,23 combining both modular hierarchical structure and a small-world property of the network (i.e. Hausdorff dimension d H = ∞). Interestingly, Complex Network Manifolds tessellate D dimensional spaces and the dimension D, together with their modular structure, strongly affects their synchronization properties. As a result of their complex network geometry, these network structures can sustain frustrated synchronization for a large range of their parameter values, with the stability of this phase depending crucially on the dimension D. The observed spatio-temporal fluctuations of the order parameter appear when the dynamics is strongly affected by localized eigenvectors, identifying "rare" regions of activity, and driving global oscillations of the synchronization global order parameter.

Complex Network Manifolds
Definition and basic structural properties. Simplicial complexes are natural objects to be considered when investigating network geometry. In fact, intuitively they can be interpreted as geometrical network structures built by geometrical "LEGO blocks" (the simplices).
A d-dimensional simplicial complex is formed by d-dimensional simplices (fully connected networks of d + 1 nodes) such as nodes (d = 0), links (d = 1), triangles (d = 2), tetrahedra (d = 3) etc., glued along their faces. Here by a face of a d-dimensional simplex, we indicate a δ-dimensional simplex with δ < d formed by a subset of its nodes. A simplicial complex has the additional property that if a simplex α belongs to the simplicial complex  (i.e. α ∈ ) then also all its faces α′ ⊂ α belong to the simplicial complex  (i.e.  α′ ∈ ). Complex Network Manifolds (CNM) 19,20,22 are generated by a non-equilibrium growing network dynamics. Specifically, Complex Network Manifolds are growing simplicial complexes of dimension d. CNM are discrete manifolds generated by gluing subsequently d-dimensional simplices along their (d − 1) faces. Every (d − 1) face α of the CNM is characterized by an incidence number n α indicating the number of d-dimensional simplices incident to it minus one. Initially (at time t = 1), the CNM is formed by a single d-dimensional simplex. At any subsequent step (at time t > 1) a new d-dimensional simplex is glued to a (d − 1) -face α with probability In ref. 22 the exact degree distribution of CNMss has been analytically derived. Mainly, the degree distribution is an exponential degree distribution for dimension d = 2 and a power-law for dimension d > 2 with power-law exponent γ given by The resulting network is small-world and has an infinite Hausdorff dimension (see Supplementary  Information). For instance, if we consider a CNM of dimension d = 3, we can embed the manifold in a d = 3 dimensional sphere of radius one without any crossing of the faces. To this end, the first tetrahedron is a regular tetrahedron inscribed into the 3-dimension sphere of radius one, and every new tetrahedron has the new node on the surface of the sphere and shares a triangular face with the already existing simplicial complex. Since in this embedding all the nodes of the CNM are on the D = 2 dimensional surface of the sphere, the network can be also interpreted as a tessellation of the D = 2 dimensional manifold formed by the surface of the d = 3 sphere (see Fig. 1).
In this D = 2 triangulated space, the elementary move of adding a new tetrahedron corresponds to the placement of a new node in a middle of a randomly selected triangle of the D = 2 triangulation and the establishment of three new links between the new node and each of the nodes of the selected triangle. Most notably, this construction reveals that CNMs in d = 3 are random Apollonian Networks [34][35][36] .
Similarly, a CNM of dimension d = 4 can be interpreted either as a 4-dimensional manifold with boundary or as the tessellation of a D = 3 space. In this case, the addition of a new 4-d simplex corresponds to the selection of a tetrahedron forming the tessellation of the D = 3 space, the placement of a new node in the middle of it and the establishment of four new links between the new nodes and each of the nodes of the selected tetrahedron.
Therefore, we can naturally associate both dimension d and dimension D of its natural embedding spaces to the CNM formed by simplices of dimension d.
A further dimension of any CNM is its spectral dimension characterizing the property of its Laplacian spectrum. This dimension will be characterized in the next paragraph.
Spectral and localization properties. The spectral properties of complex networks have been shown to be particularly relevant to reveal the interplay between network structure and synchronization dynamics on the network. Here, we emphasize the relationship between spectral properties and network geometry. To this end, we characterize the spectral properties of the normalized Laplacian matrix L characterizing any given CNM. The normalized Laplacian L has elements where a ij define the adjacency matrix (with 0 entries for non-connected nodes and 1 entries for connected ones) and k i the degree of node i in the associated network skeleton. The normalized Laplacian has real eigenvalues A large number of complex networks are characterized by a spectral gap, i.e. their second smallest (so-called Fiedler) eigenvalue λ 2 does not approach zero as the network size grows. On the other hand, CNMs, like regular lattices, have a spectral gap that approaches zero for large network sizes (see Fig. 2). This ensures that one can define the spectral dimension d S , which characterizes the power-law scaling of the density of eigenvalues ρ(λ) for λ ≪ 1 as 37,38 , S d s can be calculated starting from the associated cumulative distribution ρ c (λ) determining the probability that a random eigenvalue is less than λ, i.e. Remarkably, CNMs formed by d-dimensional simplices have spectral dimension    Fig. 2(a), we show the cumulative eigenvalue distribution ρ c (λ) and the fitted power-law behavior for small λ, which yields d S = 2.00 (2) The eigenvectors of the normalized Laplacian are also useful to reveal relevant properties of the CNMs. Given that the normalized Laplacian is an asymmetric matrix, it is characterized by a set of left u L λ and a set of right eigenvectors λ u R which, in general, do not coincide. Left and right eigenvectors associated with the same eigenvalue are normalized according to the condition The localization of any given eigenvector can be quantified using the participation ratio Y, which is an indicator of the number of nodes on which such an eigenvector has a significantly different from zero value, i.e. Fig. 2(b), we show that a very large fraction of eigenvectors are localized on a small fraction of nodes, as indicated by a small value of their participation ratio Y compared with the total number of nodes of the network, N. This differs from what happens in regular lattices, where eigenfunctions (oscillation modes) are typically delocalized.

Synchronization.
Here we investigate the synchronization properties of coupled oscillators -defined as in the Kuramoto model 25,26,39-41 -running on top of CNMs as a function of the network dimensionality. To each network node i, we associate a non-linear oscillator whose phase θ i changes in time according to its own internal frequency ω i and the coupling to other connected oscillators, obeying where σ is the control parameter tuning the overall coupling strength, and the internal frequencies ω i are random variables independently drawn from a normal distribution with mean 0 and variance 1, i.e. (0, 1)  . Previous studies on frustrated synchronization on complex networks have shown that whereas the particular choice of the distribution of ω i can affect quantitatively the phase transition between the asynchronous and a synchronous phases, its qualitative properties are quite robust 31 , therefore we restrict the study to this exemplary case. In order to attenuate the effect of possible heterogeneities in the degree distribution -and to emphasize that of network geometry -in Eq. (9) the effective coupling between node i and its neighbours has been normalized by its degree k i . In order to quantify the degree of synchronization, we employ the standard Kuramoto order parameter, defined as where R(t) ∈ [0, 1] is a real variable (function of t) that quantifies the level of global synchronization, and φ(t) gives the average global phase of collective oscillations 2,25,39 .
Several previous works have analyzed the effect that the underlying network topology has on the synchronization properties of Kuramoto oscillators. For instance, in fully connected networks, as well as in random Poissonian networks, the dynamics of the Kuramoto model yields a continuous phase transition from an incoherent state (with  R 0) to a coherent one (with  R 1) at the critical coupling σ c 39-41 . In regular lattices of dimension d it has been shown that global synchronization is only possible for d > 4; in dimensions 2 < d ≤ 4 only entrained frequency synchronization, but not phase synchronization, is observed, whereas in dimension d ≤ 2 synchronization is not observed 42,43 . Furthermore, as mentioned above, it has been recently shown that also a regime of frustrated synchronization -akin to Griffiths phases 30 -characterized by global oscillations of R(t), can emerge in complex networks with hierarchical and modular structure 31 . This regime has been previously associated only with large-world networks, i.e. networks with a finite Hausdorff dimension 30 .
In what follows we illustrate that it is the network spectral dimension what mostly determines the resulting synchronization properties, so that small-world networks -with an infinite Haussdorf dimension -of which CNMs are an example here, can also display a regime of frustrated synchronization.

Synchronization and Frustrated Synchronization.
We have performed numerical analysis of the Kuramoto dynamics on CNMs which reveals, for a wide range of coupling values, a frustrated synchronization phase in which the global order parameter has large temporal fluctuations. In Fig. 3 we show the global order parameter R(T) calculated at large times T for a given fixed CNM of size N and given internal frequencies In order to characterize the frustrated synchronization phase and to assess whether a true synchronization transition is observed for CNMs of D = 2 and D = 3, we have performed an extensive computational analysis of the considered Kuramoto dynamics averaged over different CNMs and different realizations of the internal frequencies. As a way to characterize network oscillations, in Fig. 4 the average order parameter R and its standard deviation σ R , calculated after a transient time, are shown as a function of the coupling constant σ. The large values of the standard deviation indicate the region of the phase space in which frustrated synchronization is  observed. Our finite size analysis (see Fig. 4) reveals the strong influence of the CNM dimension on the macroscopic dynamics. In D = 1 global synchronization is never achieved for large network sizes, indicating that in the large-size (thermodynamic) limit synchronization is impossible. On the other hand, synchronization in CNMs with D = 2 and D = 3 is possible for small networks but gets delayed to higher couplings for increasing system sizes, and a much broader regime of large fluctuations is observed in the D = 2 case.
These numerical observations can be understood in connection with the spectral dimension of the corresponding CNM. In fact, by linearizing the Kuramoto dynamics, it is possible to extend the results obtained in 42 for regular lattices to complex networks with finite spectral dimension. These theoretical considerations reveal that for d S ≤ 2 networks cannot synchronize, whereas for d S > 4 there is always a critical value of the coupling above which synchronization is possible. Finally, for 2 < d S ≤ 4 global synchronization is not possible but an entrained synchronized state can be observed.
Interestingly, CNMs show that it is possible to realize tessellations of a D = 3 space that have spectral dimension d S = 4. These networks have the critical spatial dimension for the onset of the synchronized state. This suggests that it could be possible to have marginal synchronization in three-dimensionally embedded networks.

Spatio-temporal fluctuations of the order parameter.
Here, we investigate further the regime of frustrated synchronization by characterizing the spatio-temporal fluctuations of the order parameter. As a matter of fact, CNMs are characterized by significantly modular structures that can be revealed by community detection methods (see Supplementary Material). These communities define regions of the networks that include a set of nodes close in the embedding space. Additionally, these nodes are also more densely connected with each other than with other nodes of the network. In order to give a visual representation of these communities, in Fig. 5 we visualize single instances of CNMs in D = 1, D = 2, D = 3 and plot, by employing different colors, the communities found by using a standard (Gen-Louvain) algorithm for community detection 44,45 .
It is, thus, natural to explore the differences between the dynamical state of these communities. For that, we define a "mesoscopic" synchronization parameter as where  n is the set of nodes in the n th community and  n | | indicates the total number of nodes in the community. Here, mod n mod n , , is a real variable taking values in the range [0, 1]. Figure 6 displays the trajectory of Z mod (t) in the complex plane for some modules of exemplary cases of CNMs in D = 1, 2, and D = 3 for coupling values that maximize σ R as indicated in the caption. In these plots, a circular trajectory describes a situation in which R mod (t) is constant in time, with full synchronization of the module corresponding to R mod (t) = 1. Random trajectories around (0, 0) describe unsynchronized modules. Partially synchronized modules, on the other hand, may describe more complex, i.e. chaotic, trajectories. We can distinguish between trajectories that oscillate within a circular crown of relatively large radius, in which R mod (t) oscillates between different states of partial synchronization, and trajectories that visit the center (unsynchronized state) too. Figure 6 also displays the time series R mod (t) and the spectral decomposition S(f) of the temporal series R mod (t) − 〈R mod (t)〉 for the considered exemplary modules. Interestingly, the spectral decomposition shows that, whereas some frequencies are in fact dominant, different modules can oscillate in quite different ways, indicating diverse synchronization states. This numerical analysis reveals the relevant spatio-temporal fluctuations observed in the frustrated synchronization phase where different modules synchronize at different frequencies. Therefore, CNMs are stylized network geometries that might provide an important framework to characterize the geometrical network properties of neuronal networks and their role in brain dynamics.
Communities and localized eigenvectors. In order to reveal the relation between the community structure of CNMs and their spectral decomposition, we have characterized the localization properties of the eigenvector corresponding to the eigenvalue λ on different network communities. To this end we have evaluated for each eigenvector λ the community participation ratio Y Q , defined as where  n indicates the n th community and C is the total number of communities. The community participation ratio Y Q ∈ [1, C] indicates the number of communities on which the eigenmode is localized. Figure 7 indicates that a large number of eigenvectors are localized in one or few communities. Thus, diverse modes are activated at different communities (and at different local coupling strengths), justifying the emergence of local patches of synchronization and overall frustrated synchronization. This therefore explains the fact that in the frustrated synchronization phase we observe a dynamics highly correlated with the modular structure of the network. Coarse graining of the frustrated synchronization dynamics. In brain research brain activity is typically measured by coarse graining the neural population dynamics at the level of macroscopic brain regions. The main starting point of such research is usually the extraction of a correlation matrix between the activity of different brain regions. Here, we investigate the role of coarse graining the frustrated synchronization dynamics on CNMs by using the local order parameter R mod,n , which describes the internal synchronization of each module of the CNM. We run twice the Gen-Louvain algorithm in order to obtain a fine partition into modules of the CNM (with D = 2), whereas retaining high modularity. We then measure the Pearson correlation between every pair of different modules (see Fig. 8(a)) to obtain a correlation matrix. Subsequently, we coarse grain further these modules (reducing the modules to about one half) by considering the aggregation generated by running a single linkage clustering to the correlation matrix (see Fig. 8(b)). We observe that the modularity of the coarse grained partition remains high, thus revealing once more the coupling between the synchronization dynamics and the network topology. Moreover, we find that the coarse grained dynamics remains non-trivial, as revealed by a complex correlation matrix, indicating a sort of invariance under network hierarchical coarse graining. A similar behavior of the dynamics under coarse graining is observed in the frustrated synchronization phase of CNMs of other dimensions D.

Conclusions and Discusion
In conclusion, this study shows the rich interplay between network geometry and synchronization dynamics by investigating the Kuramoto model running on top of Complex Network Manifolds. These networks define discrete manifolds in dimension D with the small-world property and highly modular structure, and provide an ideal theoretical setting to explore the interplay between network geometry and brain dynamics. When the  The dynamical state of each community n is measured using R mod,n . The fine partition formed by 79 communities is generated by running twice the GenLouvain algorithm, yielding a modularity Q = 0.66. The coarse grained partition has been extracted by running a single linkage clustering on the correlation matrix and cutting the resulting dendrogram in order to get 30 communities. The resulting modularity of the coarse grained partition is Q = 0.60.
SCIEnTIFIC REPORTS | (2018) 8:9910 | DOI:10.1038/s41598-018-28236-w synchronization of CNMs of dimensions D = 3 and D = 2 is compared, it is observed that both network structures sustain a regime of frustrated synchronization where spatio-temporal fluctuations of the order parameter are observed. In this regime, network modules are characterized by different synchronization frequencies. However, finite CNMs in D = 3 are much more favourable to sustain synchronized states than CNMs in D = 2. These results help shedding light on the experimental finding that 3D scaffolds favor neuronal network dynamics, as recorded in calcium activity experiments, with respect to 2D geometries for neuronal cultures, and that relevant features of brain dynamics are a consequence of its 3D topology, as experimentally observed in ref. 17 . Moreover, our study shows evidence that CNMs that can be embedded in dimension D = 3 may have spectral dimension d S = 4, i.e. the critical spectral dimension for the onset of a global synchronous phase, allowing us to observe both a fully synchronized regime and a frustrated synchronization regime on a finite network. Also, our work reveals that non-trivial synchronization states can emerge even in small-world networks, with an infinite topological dimension.
On a wider perspective this work reveals the important role of the spectral dimension and the localization of eigenvalues in promoting the frustrated synchronization phase and opens new research lines to relate network geometry and brain dynamics.

Methods
Computational analyses are performed in MATLAB software. Integration of system dynamics are carried out using the ode45 function, which uses a non-stiff 4-th order integration algorithm with adaptive steps. Modularity analyses are performed via the Generalized Louvain algorithm 44,45 . The codes for generating the Complex Network Manifolds, which are equivalent to the Network Geometry with Flavor s = −1 19 , can be found in this public repository 46 .