Pattern invariance for reaction-diffusion systems on complex networks

Given a reaction-diffusion system interacting via a complex network, we propose two different techniques to modify the network topology while preserving its dynamical behaviour. In the region of parameters where the homogeneous solution gets spontaneously destabilized, perturbations grow along the unstable directions made available across the networks of connections, yielding irregular spatio-temporal patterns. We exploit the spectral properties of the Laplacian operator associated to the graph in order to modify its topology, while preserving the unstable manifold of the underlying equilibrium. The new network is isodynamic to the former, meaning that it reproduces the dynamical response (pattern) to a perturbation, as displayed by the original system. The first method acts directly on the eigenmodes, thus resulting in a general redistribution of link weights which, in some cases, can completely change the structure of the original network. The second method uses localization properties of the eigenvectors to identify and randomize a subnetwork that is mostly embedded only into the stable manifold. We test both techniques on different network topologies using the Ginzburg-Landau system as a reference model. Whereas the correlation between patterns on isodynamic networks generated via the first recipe is larger, the second method allows for a finer control at the level of single nodes. This work opens up a new perspective on the multiple possibilities for identifying the family of discrete supports that instigate equivalent dynamical responses on a multispecies reaction-diffusion system.

Networks of interacting elements provide a formal representation to model the structure of many complex systems 1,2 , such as social interactions 3 , transportation models 4 , ecological systems 5 , and neuronal functions 6 . The dynamical behavior of such systems cannot be explained as a simple superposition of the dynamics of the single units: one also needs to account for how the different elements interact (type of coupling) and which connectivity structure links the elements (network topology). Dynamical systems on complex networks have been widely studied 7 . However, a general theory to fully resolve the subtle interplay between network structure and ensuing dynamics is still lacking, notwithstanding notable efforts which combine fundamental [8][9][10] and applied expertise 11,12 . In the particular case of reaction-diffusion systems with identical single node dynamics, the coupling between the units always admits a homogeneous state where all the nodes are synchronized. If such equilibrium proves unstable, any arbitrary perturbation grows through the unstable directions giving rise to irregular spatio-temporal patterns [13][14][15][16] . This irregular behaviour represents the functional response to a given input which can be traced back to the structural characteristics of the underlying network 17 .
Whereas, in general, a random change on the network topology would lead to a modification of its dynamical response, it is also presumable that a specific outcome is not unique from a particular network, so that similar patterns might arise from different isodynamic topologies. The identification of different compatible structures that give rise to the same dynamical behavior represents an important leap forward in the study of complex networks dynamics. For instance, it paves new roads to devise network reconstruction protocols, in cases where more than one network can correspond to the same dynamical output. Upon analyzing the common topological properties that contain the dynamical information of the system, one might overlook the inaccessible details of the topology. In this work we propose two different methods to, given a specific network, generate a second network isodynamic to the first one, so that they share the same dynamical response. A method to generate isodynamic networks for the heterogeneous Kuramoto model was already introduced in 18 . Here we assume instead a generic class of systems characterized by identical single-node dynamics and nonlinear diffusive coupling. The system admits a homogeneous equilibrium solution, whose stability reflects the underlying network topology as stored in the spectrum of the associated Laplacian operator 19,20 . We exploit the idea that, in the vicinity of an unstable (homogeneous) equilibrium, only a limited subportion of the overall network architecture proves significant for the emergence of a specific dynamical behavior [21][22][23] .
The first method is exact, since it directly acts on the subspace generated by the Laplacian eigenvectors characterizing the stability of the equilibrium state. On the other hand, the second technique relies on a Monte-Carlo algorithm that allows for a neater control in terms of network topology. In brief, the goal of this work is to prove that it is possible to provide two (or more) structurally different complex networks so that inserting the same input on them can lead to the same output.

Dispersion Relation and Pattern Formation
Let us consider a generic system composed of N identical entities linked through a complex network. At time t the activity of node j is described by an m-dimensional variable w j (t) ∈  m . Starting from a specific initial state w j (0), the dynamics of w j evolves according to and    → : m m are generic continuous functions. The first term on the right hand side of the equation describes the self-dynamics of each individual node j, and it is denoted as reaction term. The second term instead, defines the interaction of node j with the other nodes of the network: the existence of a coupling is set by the adjacency matrix A, while the interaction shape is established by function . Our analysis requires two assumptions: (i) there exists an equilibrium w j = w * for the uncoupled problem = w w ( ) j j  , and (ii) function  annihilates in zero ( = 0 0 ( ) ). The uncoupled equilibrium w * can be either a stationary fixed point or a limit cycle. Whatever the case, the second condition ensures that w * becomes a homogeneous solution of the system (1). A wide class of models corresponding to this setup are reaction-diffusion systems, where  is a linear function. This constraints can be also generalized in order to include networks of homogeneous Kuramoto-Daido systems with arbitrary coupling functions [24][25][26] .
At the homogeneous equilibrium w * no information is flowing through the network. However, if such state proves unstable, a generic small perturbation δw of w * can develop to an irregular spatio-temporal pattern. A linearization of equation (1) around w * provides the time evolution of the perturbation δw j ∈  Nm as a system of Nm linear ordinary differential equations, where J ∈  Nm ×  Nm is the Jacobian matrix of the system. Therefore, the stability analysis of w * turns into studying the high-dimensional operator J(w * ). Nevertheless, it is possible to link the diagonalization of the Jacobian to that of a simpler operator. In the linearized regime, the flow of the quantity δw j to another node i of the underlying network is described by is the in-degree of node i and δ ij the Kronecker delta. The resulting N × N matrix Δ is known as Laplacian of the network. The Laplacian matrix always has a zero eigenvalue associated to a uniform eigenvector φ , and all remaining eigenvalues are negative.
By expressing δw j in the basis of the Laplacian eigenvectors Φ = (φ (α) ), and making use of the corresponding eigenvalues Λ = (Λ (α) ) it is possible to decouple the Nm equations from system (2), thus reducing the problem to N uncoupled m-dimensional systems indexed by α (see Supplementary information). Each of the reduced systems is described by the reduced Jacobian If J α is time independent, as when w * is a fixed point, the stability analysis is simply assessed by its m eigenvalues λ λ =  [27][28][29] , which represent the analogue of λ (α) for a time-dependent Jacobian (see Supplementary information). In this case, the exponents controlling the stability of w * are given by λ(α) : = max k µ α ( ) k ( ) In both instances, λ provides the stability of w * as a function of the Laplacian eigenvalues. Such relation takes the name of dispersion relation. If the dispersion relation of a subportion N c of the total N eigenvalues Λ (α) is positive, then any perturbation would grow through the unstable modes. The resulting irregular spatio-temporal pattern thus represents the unpredictable response of the system to a specific input signal. It is however reasonable to suppose that most of the relevant information is stored into the unstable manifold of the homogeneous solution w * , as we will prove later on.

Topology Modification
In this section we propose two methods to modify a given discrete support of the system while preserving the relevant directions for the emergence of the pattern. suppose, without loss of generality, that the subset of modes to be left invariant are the first n. We define the diagonal matrix Λ ∼ , such that The corresponding eigenvectors are modified by performing a change of basis where I n is the identity matrix of dimension n and R N−n ∈SO(N − n) is a random rotation matrix. In order to practically obtain such rotation we perform QR decomposition of a random matrix. Taking care of preserving the uniform mode, which corresponds to an identical perturbation acting independently on each node, the transformation Local rewiring. The previous method does not provide any control whatsoever on the topological modifications introduced in the new Laplacian. For this reason, we propose a secondary route that acts at the level of single nodes. For many network structures, the Laplacian eigenvectors are well localized on the network, i.e., their coordinates in the original vectorial space mostly involve a small subset of nodes, different for each eigenvector 30 . Therefore, modifying some eigenvectors means acting mainly on the connections of a specific subnetwork, and vice versa. With this method we aim to identify and modify links among nodes that are poorly involved on the n-dimensional manifold we want to preserve.
We rely on a Monte-Carlo algorithm that proceeds as follows. Given the original network, we choose a random non-diagonal entry of the adjacency matrix, A ij . If the entry indicates the existence of a link between node i and node j, such link is removed, otherwise it is created. We then compare the n eigenvalues and eigenvectors of the modified network with those of the original one, and only if they prove similar according to a chosen threshold τ the change is accepted (see Supplementary information for further details). The process is repeated selecting new random entries over the modified adjacency matrix until a desired number of links have been changed or the method fails to detect new entries that lead to small error.

Results and Discussion
The Ginzburg-Landau model. In order to test the different methods we use the Stuart-Landau oscillator as a reference system for the dynamics of each node. Such model is the canonical form of the Hopf bifurcation, i.e., it describes a generic limit-cycle at the vicinity of a bifurcation where a fixed point becomes unstable through a pair of complex conjugate eigenvalues crossing the imaginary axis. Therefore, it has been studied in a wide number of applications where periodic oscillatory behavior is relevant.
We consider an ensemble of Stuart-Landau oscillators occupying the nodes of an undirected network with a diffusive coupling (for a case with a directed network see Supplementary information). The equation governing the dynamics of each node j, which takes the name of complex Ginzburg-Landau (CGL) equation [31][32][33] , is where c 1 ,c 2 ,K∈ and w j ∈. Such system always accepts a homogeneous time-dependent solution, corresponding to the limit-cycle of a single decoupled element of the network We focus our analysis on this globally synchronous state.
Even though w * is time-dependent, the Jacobian J α is constant in time, and following the procedure indicated in the Supplementary information, we obtain the dispersion relation characterizing the linear stability analysis of the limit cycle solution, Depending on the system parameters, function λ might have a positive part, corresponding to the smaller Laplacian eigenvalues in absolute value. Therefore, a particular network with N nodes might have N c < N Laplacian eigenvalues associated to unstable modes. From now on, we assume that the Laplacian eigenvalues are sorted in descending order, Λ (0) = 0 > Λ (1) > … > Λ (N) , so that the N c unstable modes correspond to the Laplacian eigenvalues with index ranging from 1 to N c .
Eigenmode randomization. Erdös-Rènyi networks. We first test the eigenmode randomization method using as original topology an Erdös-Rènyi (ER) network with N = 100 nodes and average node degree 〈k〉 = 3.5. For system parameters K = 1, c 1 = 1, and c 2 = −3, the dispersion relation associated to this network has N c = 44 modes corresponding to unstable directions (see blue open circles in Fig. 1). We integrate the system using as initial condition a randomly generated small perturbation of the homogeneous limit-cycle solution. After a short transient, the system reaches a dynamical regime characterized by an irregular spatio-temporal pattern depicted Scientific RepoRts | (2018) 8:16226 | DOI:10.1038/s41598-018-34372-0 in Fig. 1(g). Notice that throughout the paper we focus only on each oscillator modulus, since the frequency of rotation does not seem to contain relevant spatial structure (see Supplementary information for more details).
We aim to generate a new topology that reproduces such dynamical behavior by leaving invariant a subset of the original network modes. Taking into account that the Laplacian eigenvalues are sorted in descending order, we preserve the first n modes and modify the rest using the eigenmode randomization method. As an illustrative example, in Fig. 1(a) we show the pattern resulting from a network where only the first n = 10 modes of the original topology are preserved and all the others have been randomized. Manifestly, such outcome has little similarity with that of the original network (cf. with Fig. 1(g)). In Fig. 1(b) we plot a zoom on the unstable part of the dispersion relation for the modified network (red dots), and that of the original topology (blue open circles). Overall is plausible to explain the disagreement between the dynamical behavior of the two networks from the large difference between the corresponding tangent space of the synchronized solution. Repeating the procedure with n = 40 preserved modes instead, leads to the pattern from Fig. 1(d). In this case, the behavior of the modified topology resembles much better that of the original network, an observation which is to be traced to the similarity between the dispersion relations of the respective unstable manifolds (see Fig. 1(e)). In order to make the comparison between different patterns more transparent, we compute the time-average modulus of each node, 〈|w j |〉 after discarding a transient of 1000 time units. In Fig. 1(c,f) we report the resulting mean node activity of the modified networks versus that of the original one for n = 10 and n = 40 respectively. An outcome from two identical patterns would lie exactly on the diagonal, while two completely independent processes would provide a random collection of points. One can then quantify the similarity between the modified and original patterns in terms of squared Pearson correlation coefficient R 2 , which is 0.38 for n = 10 and 0.97 for n = 40. In Fig. 1(h) we show the outcome of repeating this analysis systematically for increasing number of preserved modes n. When the modes to be preserved are selected to be the first n, the correlation between the patterns of the modified and original networks increases quickly with n, reaching a very good agreement when approximately all the 44 unstable modes are preserved (see red circles). On the other hand, if the n invariant modes are randomly selected among all the N, even when a large number of directions are maintained, there is no guarantee that the modified network will respond similarly to the original one (see gray triangles and squares).
The eigenmode randomization technique does not provide, a priori, any information about the structure of the resulting network. In Fig. 2(a,b) we show the topologies of the original network and a modified version where all the 44 unstable modes have been preserved. The first obvious difference is the existence of weighted links, including negative ones (see blue edges), which are absent in the original network. The adjacency matrix loses its sparsity and the network becomes highly connected, although most of the new links are very weak. In order to compare the two networks we focus on the strength s j of each node, defined as = ∑ = s A j m N jm 1 , which extends the concept of node degree to weighted topologies. In Fig. 2(c,d) we plot the degree distribution of the original network, and the strength distribution of the modified topology, respectively. Although the distribution of the new network looks different from the original one, it does not present any particular structure, as expected for a random topology. In fact, the average strength of the new network coincides with the average degree of the initial support.
Scale-Free networks. In order to further investigate the relation between the Laplacian eigenmodes and network topology, we next move to a case where the network has a specific structure. For this purpose we use Scale-Free networks (SF), characterized by a typical power law degree distribution. In Fig. 3(a) we show an example of such degree distribution for a network composed of 150 nodes built using the Barabasi-Albert 34 preferential attachment algorithm. Fixing the parameters of the CGL equation as K = 1, c 1 = 1.2, and c 2 = −10, the dispersion relation of this network shows N c = 76 unstable modes. We apply the eigenmode randomization technique to generate a new support taking care to preserve all the unstable directions, but none of the stable (n = N c ). The patterns resulting from inserting the same perturbation to the synchronized solution on both networks are highly alike, as can be seen from the correlation Fig. 3(d). Differently from the ER case, here the modulus of each node gets stationary after some transient, so such time-average activity is free from statistical fluctuations.
The topology of the new network is, again, highly connected and involves negative links. Morevoer, the node strength distribution (see Fig. 3(b)) does not preserve the scale-free structure of the original topology. To identify the nature of this strong changes, in Fig. 3(c) we plot a one to one comparison between the degree of each node of the original support and the corresponding strength on the new topology. This analysis reveals that nodes with smaller degree on the original setup have a comparable strength after the modification, whereas the hubs become much weaker.
The localization properties of the Laplacian eigenvectors allow us to understand this situation. Inspired by the analysis performed by Hata and Nakao in 30 , we depict in Fig. 3(f)   This implies that all the networks generated with the eigenmode randomization procedure will mostly differ from the original one for what concerns the nodes characterized by a large degree k j . For this reason the characteristic shape of the degree distribution is not preserved in the second network. Nevertheless, in dynamical systems with a different shape of the dispersion relation it might be possible to maintain the tail of the distribution if the stable modes are to be found among the first Laplacian eigenmodes (see a specific example in the Supplementary information).
Finally, we repeat the procedure of systematically increasing the number of invariant modes n from 1 to 150 and tracking the resulting correlation coefficient, Fig. 3(f). For each new network we analyze the dynamical patterns resulting from inserting three different randomly generated perturbations. As for the ER case, the correlation between the dynamics of the original and modified topologies increases quickly with n, being nearly optimal when all the unstable modes are preserved (see vertical dashed line). Nevertheless, the choice of the initial condition is clearly relevant: whereas in one case preserving around 40 modes already provides a good agreement (see black crosses), in other situations one might need to preserve also a considerable number of stable directions to reproduce the original pattern (see red circles).
Local rewiring. The second method we present offers the advantage of a larger control on the resulting network topology so that, for instance, one can keep the network binary. On the other hand, one needs to allow for some variability also on the unstable manifold due to the non perfect localization of the Laplacian eigenmodes on a specific subnetwork.
We first apply the local rewiring technique on a ER network with N = 400 nodes and average degree 〈k〉 = 20. For system parameters K = 0.15, c 1 = 1, and c 2 = −3, the system displays 64 unstable modes. The algorithm stops after 100 positions on the adjacency matrix have been changed. For comparison purposes, we also created three different networks where 100 entries on the adjacency matrix have been changed totally at random. Blue circles in Fig. 4(a) show the correlation between the patterns of the original and the resulting network, whereas the gray symbols correspond to that of the random rewired graphs. The dynamical behavior of the network obtained using the local rewiring method provides a better agreement with the original pattern. We obtain similar results when using SF network generated with the Barabási-Albert algorithm (see Fig. 4(b)). Although these results change upon inserting different initial perturbations, the rewiring technique outperform the random modified networks in most of the cases (see Supplementary information for additional results). In terms of topology changes, in Fig. 4(c,d) we plot the total number of links that have been modified for each node of the network. It is clear that most of the changes correspond to nodes with larger degree which, as shown in the previous section, mostly involve the Laplacian eigenvectors associated to stable directions in the case of the CGL equation.

Conclusions
In the present work we have devised two specific strategies of network generation so as to mimic the pattern obtained by a former sample network when both are subject to the same non-linear reaction-diffusion system. The first method is analytical, and provides isodynamic networks at the cost of changing important topological features of the original graph. The second technique instead makes use of a Monte-Carlo algorithm, allowing for a larger control in terms of network topology, but providing less accurate results. Both methods rely on a preliminary identification of a manifold generated by the Laplacian eigenvectors associated to homogeneous solution instability, which has to be preserved during the modifications.
This work sheds light on the strong dependence between network topology and dynamical patterns. The localization properties of the Laplacian eigenvectors unveil the existence of underlying topological features that lead the dynamical response of the network. In particular, it is possible, in some cases, to identify a subnetwork or a set of nodes that can be acknowledged as practically irrelevant for pattern formation. It is then clear that distinct networks with differences solely localized on these particular substructures can, under apt conditions, generate the same irregular patterns. As a remarkable example of this situation, we have shown that the SF networks, often taken as reference model to describe many real systems, can be, in the case of the Ginzburg-Landau system, replaced by graphs characterized by highly different structures.
Arguing on this line, the two methods proposed here may be placed as novel schemes for network generation, with a purely functional angle: instead of building the specific network structure brick by brick, paying attention to local characteristics of nodes, the complex graph is globally established and constrained to yield an a priori chosen dynamical output.
As already mentioned, the issue of generating isodynamics networks has been tackled in 18 with reference to the heterogeneous Kuramoto model. In this latter work, the focus was placed on preserving global collective Figure 4. Results of the local rewiring method applied to a ER with average degree 〈k〉 = 20 and 64 unstable modes (a,c), and to a SF network with 63 unstable modes (b,d). Both topologies consist of N = 400 nodes. At each step the tolerance error of the algorithm is τ = 0.1. In total, the new networks present 100 modified links with respect to the original topologies. (a,b) Correlation between patterns of the original and modified networks. Colored circles correspond to the results obtained using the local rewiring algorithm. Gray pluses, crosses, and stars correspond to three different networks generated by random rewiring as many links as the network outcoming from the local rewiring procedure. The squared coefficient correlations for the modified ER network is 0.98 (blue circles) whereas the randomly rewired networks provide 0.67, 0.62, and 0.94. The modified SF network has R 2 = 0.86 (red circles), whereas the random networks provide 0.49, 0.50, and 0.67. (c,d) Number of link changed for each node of the network. The nodes are sorted in descending degree order.
Scientific RepoRts | (2018) 8:16226 | DOI:10.1038/s41598-018-34372-0 variables, as e.g. reflecting the displayed level of synchronization. At variance, we here deal with a local definition of "dynamical invariance", which extends to individual nodes. In carrying out the analysis, we have however postulated the existence of an underlying homogeneous equilibrium, a working hypothesis which makes the method unsuited to deal with quenched disordered, as in 18 . Extending the proposed technique to include the case of heterogeneous equilibria is an intriguing direction of investigation which is left for future work.