Synchronization Invariance Under Network Structural Transformations

Synchronization processes are ubiquitous despite the many connectivity patterns that complex systems can show. Usually, the emergence of synchrony is a macroscopic observable, however, the microscopic details of the system, as e.g. the underlying network of interactions, is many times partially or totally unknown. We already know that different interaction structures can give rise to a common functionality, understood as a common macroscopic observable. Building upon this fact, here, we propose network transformations that keep the collective behavior of a large system of Kuramoto oscillators functionally invariant. We derive a method based on information theory principles, that allows us to adjust the weights of the structural interactions to map random homogeneous -in degree- networks into random heterogeneous networks and vice-versa, keeping synchronization values invariant. The results of the proposed transformations reveal an interesting principle; heterogeneous networks can be mapped to homogeneous ones with local information, but the reverse process needs to exploit higher-order information. The formalism provides new analytical insight to tackle real complex scenarios when dealing with uncertainty in the measurements of the underlying connectivity structure.

The study of dynamical processes running on top of complex networks has become a central issue in many research fields, ranging from the microscopic realm of genes and neurones to the large realm of technological and social systems [1][2][3][4][5][6][7]. The interplay between topology and dynamics is crucial here to understand the physics of those complex systems under analysis. However, many times the information we can accede to about the actual topology of interactions is somehow incomplete, because of experimental limitations or because of lags on the details of the system [8][9][10][11]. Moreover, given that the only reflection of the dynamics on networks is usually a certain macroscopic observable, it turns out that many topologies are compatible with the same dynamical output, raising the problem of multi-valuation [12][13][14] (i.e. different topologies with the same dynamical response).
Following this perspective, we analyze the relation between function and structure in a novel mapping problem. Essentially, given a certain network structure and a dynamical process on top of it, we wonder how to transform the network into a different structural connectivity so that the collective behavior (i.e. the function) remains invariant. Such transformation must adjust the weights of the interactions in the new configuration to achieve the goal of having an equivalent steady-state functionality to the original structure. In this letter, we present a new formalism, based on the maximum entropy principle [15,16], to derive analytical transformations for the resulting weights when only local information (at the nodes' scale) is available. Furthermore, we show that the mapping of homogeneous networks into heterogeneous ones is usually less accurate and requires more -costly-microscopic information than the reverse process, unveiling a symmetry-unbalance phenomenon that emerges from the partial impossibility of preserving the local structural constraints.
To derive the network transformations, we focused on a particular dynamical process, the synchronization of coupled phase oscillators. This paradigmatic example of emergent phenomena has been extensively studied [4,17,18], to unveil fundamental aspects related to the mapping problem, such as the inference of structure from response dynamics [10,[19][20][21], the dependency of the collective behavior on the topology [18,22,23] and the network optimization to maximize the stability of the fully-synchronized attractor [24,25]. The Kuramoto model (KM) [26] consists of a population of N coupled phase oscillators that evolve in time according to the set of equationṡ where θ i is phase of the i-oscillator, ω i its natural frequency, drawn from a probability distribution g(ω), λ ij are the elements of the coupling matrix Λ that capture the presence of a connection and its intensity and K is a constant coupling strength that scales all the weights. The collective behavior of the KM is usually described through the complex order parameter where the modulus r measures the overall degree of synchrony and Ψ(t) the average phase of the system. Here, we assume that the macroscopic order parameter r is the only available observable from measurements, and we look for transformations of Λ that keep this observable invariant, for any value of the control parameter K.
It is well known that particular unweighted instances drawn from the same degree distribution will produce the desired invariant collective behavior [3][4][5]. We wonder if the former invariance can be achieved for weighted networks drawn from different degree distributions, preserving the number of nodes N . We consider a target network A with a given coupling matrix A, which might be non-symmetric and directed with fixed entries λ A ij , and a candidate network B, with different coupling matrix B and λ B ij . We impose transformations of B in the form where w ij are the parameters to find. Note that we can absorb the weights of B into W, keeping only the binary values b ij of the structure of B. After this simplification, the entries of the transformed network can be written as (B ′ ) ij = w ij b ij . Furthermore, we assume that the N units are distinguishable and preserve their intrinsic properties in the trans- , which ensures that we are dealing with particular instances of networks and not with averaged ensembles. Then, the condition for functional synchronization invariance can be written as where the measurements are in the steady-state, the average refers to different initial conditions, accounting for fluctuations of order 1/ √ N , and the parameters of the dynamical process ( ω, K) are fixed in both networks.
Inspired by the derivation of statistical mechanics from information theory as a particular case of statistical inference, see [16], we tackle the functional mapping defined above as an optimization problem for the unknown weights subject to structural constraints on the networks that capture our prior knowledge on the system. The key assumption here is that Eq.(3) can be achieved by imposing a local detailed balance for the main structural properties of the nodes: the overall coupling intensities received from neighbours (or input strengths [5]). For each node, we define the zero-order input strength as s and so on. For a fixed order M , the detailed balance is given by a set of N (M + 1) equations for the s (m) i . If we let q be the N-vector of ones q = (1, 1, ..., 1) ⊤ , we have are the node structural bounds in the optimization of the weights in B ′ . The local constraint (m = 0) can be written explicitly as which ensures to preserve the overall coupling in the transformation ( i j λ A ij = i j w ij b ij ). The ansatz of Eq.(5) relies on the weighted annealed approximation [27,28], that assumes statistical similarity among nodes with the same s (0) i . This description is known to be valid in the linear regimes of the diffusion of random walkers [28] and the Master Stability Function (MSF) [24,25]. Here, if the coupling strength K is sufficiently large, Eq.(1) can be linearized, and using statistical and meanfield arguments [25], the system can be uncoupled, with each unit being driven only by its input strength. The underlying assumption is that higher order constraints (m > 0) might be required when the non-linearity of Eq.(1) plays a crucial role or the connectivity patterns are highly non-trivial (heterogeneity, correlations, etc...).
We take advantage of information theory [15], to define an appropriate objective function to optimize the unknown weights. In an uncertainty scenario, the best we can do is to rely on the Maximum Entropy Principle [16]. It states that, subject to the available data (i.e. the constraints in Eq.(4)), the probability distribution which best captures our lack of information is the one that maximizes the entropy. Here, we can interpret the weights distribution in probabilistic terms, where the input strength s (0) i is the normalization condition, and we can define the entropy [15] of a node S i as a sum over the accessible states defined as those where b ij = 1, where the normalization constant has been neglected for simplicity and it is assumed that w ij ≥ 0. We can use the method of Lagrange multipliers [16] to solve this optimization problem. The lagragian function reads as where β (m) i is the m-order lagrange multiplier of inode. By optimizing Eq.(7) with respect to the unknown weights and finding the values of the multipliers, we can derive analytical expressions for the entries of B ′ . For the zero-order case (M = 0), we obtain that can be written as w This solution is very intuitive, since it homogeneously allocates the input strength of a node into the available links. The weights are therefore equal for all the incoming links of a node (w ij is independent of the node j), implying usually a non-symmetric coupling.
The solution in Eq. (8) is precisely the scheme used in [24,25] to transform a network topology into a purely homogenous one to optimize the stability of the synchronized state in the scope of the MSF. That means that the solution is valid in the linear regime, close to the synchronization attractor. However, this solution is yet to be validated in the fully non-linear regime. We simulate the dynamics of N = 2000 oscillators following Eq.(1) with fixed g(ω) ∈ (−π, π), measuring r 2 in a quasi-static process controlled by the control parameter K ∈ [0, 0.5/N ]. We propose to map pairs of uncorrelated networks drawn from different degree distributions, that range from homogeneous in degree, Erdös-Rényi networks, to powerlaw in degree networks, which are initially unweighted and symmetric. We use the model in [29], to interpolate between both degree distributions using a single parameter α. For α = 0 we have pure power-law distributions p(k) ∼ k −γ with exponent γ = 3 while for α = 1 we obtain homogeneous random networks, keeping the average degree fixed, in our case k = 10. The mapping transformation is then as follows: we fix the topologies of a network A drawn from the model for a certain value α, i.e. the target network A α , and the candidate network B α ′ drawn for another value α ′ . Then, we compute the weights, using Eq.(8), to map the candidate network into the target one and obtain the resulting T 0 (B α ′ |A α ), where the subindex of T refers to the fact that the method exploits only zero-order information.
In Fig.(1) we present the results of the transformation for the extreme cases T 0 (B 0 |A 1 ) and T 0 (B 1 |A 0 ). The results evidence that the functional invariance is attained in the linear regime (K ≫ K c ) for both transformations. However, there is a clear discrepancy in the transformation T 0 (B 1 |A 0 ), i.e. from a homogeneous in degree network towards an heterogeneous, power-law, network. This discrepancy shows that, when Eq.(8) is applied, homogeneous networks are not able to capture the role of heterogeneous connectivity patterns. To improve the accuracy of the T 0 method in the mapping, we need to include higher-order constraints. We extend the detailed balance to a further order (M = 1) by imposing that, for each node, the transformation must also preserve the first-order input strengths s Note that s (0) j is the same at both ends of Eq.(9) because we still retain the constraint presented in Eq. (5). We aim to maximize Eq.(6) subject to Eq.(5) and Eq.(9). The lagrangian in Eq. (7) can be written explicitly as By imposing dL/dw ij = 0 and isolating the unknown weight w ij , we obtain the implicit expression The values of the multipliers β i are found by substituting Eq.(11) back in Eq.(9) and numerically solving the resulting system. However, the existence of real and nonnegative solutions cannot be ensured a priori. Indeed, the structural bounds are easily estimated by considering the worst-case scenarios, i.e.
The inequality in Eq.(12) turns out to be unfeasible for most nodes if the reference network is very heterogeneous in local input strength. Let us illustrate this by considering, on one hand, that A follows a power-law distribution with p(s) = cs −γ . Then, if network B is sufficiently wellconnected (k B i ≫ 1 ∀ i ∈ N ) and assuming N large, we can approximate the constraints by The first integral can be written as the Gamma function e −x x −γ dx = Γ(1 − γ). Using the well-known property Γ(z + 1) = zΓ(z) and dividing both equations, we obtain which is negative for γ = 3, thus unveiling the structural restrictions that emerge when mapping any arbitrary network into a highly heterogeneous one. On the other hand, Eq.(8) is recovered from Eq.(11) only when s (0) i ≃ s (0) , ∀ i ∈ N , i.e when A is very homogeneous in local input strength, regardless of the topology of B.
The previous reasoning unfolds the symmetryunbalance observed in Fig.(1) and suggests that the mapping can indeed be enhanced, although it is strongly limited by the structural bounds. To provide an analytical transformation that improves the performance of Eq.(8) while still preserving w ij ≥ 0, we expand Eq.(11) to first order around its average value, i.e.
j . We insert Eq.(16) into Eq.(9) to obtain an approximate value β * i ≃ β i as The solution is finally obtained by direct substitution of Eq.(17) into Eq. (11), and we denote this transformation T 1 (B α ′ |A α ). Note that T 1 does not provide uniform weighting, but depends explicitly on the balance between input strengths and heterogeneity in each node. Now we can compare the performance of transformations T 0 and T 1 in the mapping. We define, for each transformation, the dynamical error as a measure of the total difference in the synchronization diagrams between the target and transformed networks, and we define the structural error as a measure of the total difference in the first-order local structure. In Fig.(2a) we present the synchronization diagram for the extreme cases T 1 (B 0 |A 1 ) and T 1 (B 1 |A 0 ) in the same set up as before (N = 2000). We can observe a significant improvement in the transformation T 1 (B 1 |A 0 ) with respect to the zero-order method in Fig.(1), although there still are non-vanishing errors around the critical point due to the unfeasible structural bounds of Eq.(12). In Fig.(2b), we plot the dynamical σ d and structural σ s errors for different values of the parameter α in T (B α |A 1−α ). Note how the accuracy of the transformations is enhanced by T 1 for any value of α, and it is associated to a decrease in the structural error, thus validating the main assumptions of our approach. Furthermore, the approximate solution of Eqs. (11,17) can still be improved by i) considering higher-order constraints (M > 1), but then the system would become coupled and it should be solved simultaneously for all nodes, ii) extending the expansion of Eq.(11) with additional terms, iii) allowing the presence of negative interactions or indistinguishable units (without labelling the nodes in the transformation), and also iv) imposing global constraints instead of local ones (requiring costly numerical methods and global objective functions [30]). Summarizing, we have presented an analytical methodology that successfully produces functional synchronization invariant networks for the KM, by transforming the weights of the interactions, while preserving the underlying topologies, and exploiting only local structural information. We have shown that different microscopic configurations can produce the same macroscopic dynamical observables if the weights are adjusted in a way that the main local properties of the nodes are preserved. Furthermore, we have unveiled that the mapping of homogeneous networks into heterogeneous ones requires to exploit additional (up to first-order) information and it is more complicated than the reverse process, due to intrinsic structural limitations of the networks.
The presented formalism can be applied in a wide spectra of problems beyond the mapping scenario. Our framework provides a more comprehensive understand-ing of the collective behavior of oscillators on weighted and directed networks from a local perspective and can be used to make analytical predictions on them (when transformed to unweighted structures) [18,23]. Also, the transformations can induce specific features of heterogeneous networks in homogeneous ones and vice-versa, without changing the underlying structure. Straightforward examples include the possibility to induce explosive transitions in homogeneous networks (by correlating the intrinsic frequencies with the input strengths [31]) and to control the critical point of a macroscopic phase transition [3,18] only by a local readjustment of weights. From a theoretical point of view, our results are sheltered by previous works that explore information-theoretic tools to study the structure of complex networks [32][33][34] and to tackle reconstruction problems [35][36][37]. Nevertheless, here we introduce a novel connection between purely structural constraints and collective dynamical behavior. This new connection can help in refining state-ofart inference methods with driving signals [10,11] (by inferring appropriate network candidates from the available structural and dynamical information), it deeps our understanding on findings that relate weighted, directed and inhibitory interactions to optimal synchronization performance [38][39][40], and provides a new approach for evolving networks models [3,5,18], in which a network of biological units might evolve, due to an evolutionary pressure, towards heterogeneous structures that maximize the number of accesible transformations and, consequently, their potential dynamical range [41].