Relay synchronization in multiplex networks

Relay (or remote) synchronization between two not directly connected oscillators in a network is an important feature allowing distant coordination. In this work, we report a systematic study of this phenomenon in multiplex networks, where inter-layer synchronization occurs between distant layers mediated by a relay layer that acts as a transmitter. We show that this transmission can be extended to higher order relay configurations, provided symmetry conditions are preserved. By first order perturbative analysis, we identify the dynamical and topological dependencies of relay synchronization in a multiplex. We find that the relay synchronization threshold is considerably reduced in a multiplex configuration, and that such synchronous state is mostly supported by the lower degree nodes of the outer layers, while hubs can be de-multiplexed without affecting overall coherence. Finally, we experimentally validated the analytical and numerical findings by means of a multiplex of three layers of electronic circuits.

dependencies of the phenomena are studied, using perturbation stability analysis. The robustness of the relay synchronization against de-multiplexing the layers is reported, revealing the key role of low degree nodes in maintaining the layers coordination. Finally, the findings are experimentally validated in a multiplex network of electronic circuits.

Results
Model. We start by considering 2M + 1 layers (or networks), arranged as shown in Fig , with  ∈ u i k m , i = 1, …, N, and whose intra-layer interactions are encoded through the Laplacian matrices   = { } k ij k . The layer stack is symmetric with respect to k = 0 in such a way that Laplacians  k and −k  have the same structure. The dynamical systems are also paired: nodes belonging to layers U +k and U −k are identical to each other, and different (in some parameter) from the rest of the layers. Consequently, layer k = 0 has no counterpart, and acts as a relay between all layers situated above and below it.
Layers are coupled in a multiplex configuration, and the dynamical evolution of the system is described by the following equations: (with   → f : k m m representing the vectorial functions evolving each dynamical unit), are identical for the same |k|. G, H are the m × m matrices representing respectively the linear intra-(G) and inter-(H) layer coupling schemes.  N is the N × N identity matrix, σ k is the intra-layer coupling strength within layers k and −k, and λ is the inter-layer coupling strength.
Due to the reflection symmetry of the system under study (i.e. as long as the U +k and U −k layers are identical for all k), a synchronous inter-layer evolution (with layers evolving in a pairwise synchronized fashion, i.e. where U + = U −k ) at all k without necessarily implying U k = U k′ for k ≠ k′) is a solution of Eq. (1), independently of intra-layer synchronization 23 (i.e. independently on whether the state of the systems within each layers are same color) are networks of identical oscillators with the same topology  k and intra-layer coupling σ k and whose dynamical state is described by the variable U k and U −k , respectively. The multiplex is symmetric with respect to the layer k = 0 and the nodes are coupled to their replicas in the rest of layers with an inter-layer coupling strength λ.
Considering the smallness of δ δ δ δ and expanding around the inter-layer solution up to first order, one obtains a set of M linearized vector equations for the perturbations δU k : where J denotes the Jacobian operator, δ kM is the Kronecker delta accounting for the boundary condition at k = M (as the stack end layers U ±M are only connected to the previous neighbor layer). The vector k describes the dynamical state of any of the k = 0, …, M layers at the synchronous state U k = U −k ≠ U 0 and, therefore, the whole dynamics is reduced to the dynamics of M + 1 layers. Such evolution at the node level is given by: , N is the node index, k, q = 0, 1, …, M, and g(u) = Gu and h(u) = Hu are the projections of the inter-and intra-later coupling operators to the node level. Notice that, since each paired layers k and −k is k , each layer acts therefore as a relay to the rest of the stack. The problem consists now in solving MmN linear equation (2), together with solving in parallel the (M + 1)mN nonlinear equation (3)  as a function of the system parameters actually gives the necessary conditions for the stability of the synchronous solution: whenever < MLE 0, perturbations transverse to the manifold will die out, and the multi-relay synchronous solution will be stable.
In order to monitor the synchronization error between layers, we define the inter-layer synchronization error as, where ||·|| stands for the Euclidean norm and q, k are the layers' indexes, such that E −k,k denotes the inter-layer synchronization error of mirror layers. Without lack of generality, In our numerical simulations we consider two types of topologies where layers are either (i) Erdös-Renyi 39 (ER) or (ii) scale-free (SF, generated by means of the Barabási-Albert's algorithm 40 ), in all cases with N = 500. We classify the layer stacks regarding the topology sequence of each layer. For instance, a triplex where the three layers have ER topology will be denoted as EEE, and a system where two identical SF layers are mediated by a center ER will be denoted as SES. The nodes are chaotic Rössler oscillators 41 , defined by the m = 3 state vector u = (x, y, z) whose autonomous evolution is given by k k k and the heterogeneity between the layers is introduced through the parameter a k , such that each layer develops a different chaotic dynamics. In our case study, the intraand inter-layer coupling functions are set to be = = z g u Gu ( ) (0, 0, ) T and = = y h u Hu ( ) (0, , 0) T respectively. These coupling schemes ensure that intra-layer synchronization is prevented when layers are isolated and not multiplexed (class I layers, according to the standard master stability function (MSF) classification established in ref. 2 ) whereas multiplexed nodes along the layers can synchronize for a coupling strength λ above a given threshold (class II MSF) 42 .
Layers with identical topology. With the aim of determining whether relay synchronization can be achieved in a multiplex configuration let us first consider the multiplex structure defined by Eq. Results are collected in Fig. 2, where the synchronization error between the outer layers E −1,1 is plotted versus the inter-layer coupling λ for several values of the intra-layer couplings σ 1 and σ 0 in the outer and relay layers respectively, with σ 1 = σ 0 . In all cases, there is a critical coupling λ ⁎ above which complete synchronization between layers k = 1 and k = −1 occurs, that is, E −1,1 = 0 is achieved, while the relay layer (k = 0) remains unsynchronized to any of the two outer layers (k = 1, −1) as shown in the inset where > E 0 0,1 for any parameter choice. In addition, the calculation of the corresponding MLE given by Eq. (2) (bottom panel of Fig. 2) confirms that the relay synchronous solution U −1 = U 1 reaches stability (MLE < 0) at the same critical λ * where the error between the outer layers is zero, as indicated by the vertical lines. Therefore, one can conclude that inter-layer MLE is a useful tool for reducing the system's dimensionality and use it for evaluation of the critical inter-layer coupling λ * from now on.
In order to better understand the different roles played by external and relay layers, we show in Fig. 3 the critical inter-layer coupling value in the parameter region (σ 0 , σ 1 ), that is, when the intra-layer coupling σ k is different SCIeNtIfIC REpoRtS | (2018) 8:8629 | DOI:10.1038/s41598-018-26945-w for the relay and outer layers. It can be seen that the system's ability to synchronize is practically unaltered with σ 0 , while increasing σ 1 makes the value of λ * to drop drastically. This therefore reveals that multiplex relay synchronization is much more sensitive to changes affecting the mirror layers than to those arising in the transmission layer.
Our results can be generalized to any number of layers. As an example, we report also the case M = 2, which corresponds to two outer layers above (k = 1, 2) and below (k = −1, −2) the relay layer (k = 0). We choose a −1 = a 1 = 0.2 and a −2 = a 2 = 0.3, and a 0 = 0.25 for the central layer. The results stand for any other parameter choice. In Fig. 4 we plot the inter-layer synchronization errors E −1,1 (void markers) and E −2,2 (full markers), vs. the inter-layer coupling λ for several values of the intra-layer coupling σ. As in the triplex case, the critical λ * at which complete inter-layer synchronization is achieved depends on σ, but it is the same for both pairs of layers, as E −1,1 and E −2,2 drop to zero simultaneously. In the inset we plot the inter-layer synchronization errors between the non-paired layers, E 0,1 , E 1,2 to check that they remain mutually incoherent. Therefore, we have verified that relay synchronization also occurs in cascade for arbitrarily high-order multiplex systems, provided a structural and dynamical symmetry is conserved.
Layers with non-identical topology. So far, we have dealt with multiplexes of pairwise identical layers.
However, this condition is too strong a limitation to hope that it would capture and properly represent the case of many real systems. The next step needed for generalization is studying then the relay synchronization scenario in  the case in which the topology of the relay layer differs from that of the outer layers. In Fig. 5 we have reported the critical inter-layer coupling λ * in two heterogeneous triplex cases: (a) a pair of Erdös-Rényi layers mediated by a scale-free relay layer (ESE situation) and (b) the opposite case where SF layers are connected through a ER layer (SES). Each case is compared with the topologically homogeneous EEE and SSS structures, respectively. For the sake of simplification and of better assessment of the role of the topology, we keep σ 0 = σ 1 . Figure 5(a) shows that, for a large range of intra-layer couplings, the mediation of a SF relay facilitates the synchronization between the paired layers, since λ * in the ESE case (void blue circles) is smaller than the one corresponding to the homogeneous case (EEE, full blue circles). On the contrary, a relay ER layer intermediating between two outer SF layers ( Fig. 5(b)) does not determine a significant difference as long as the intra-layer coupling strength is low, but increases the threshold λ * for higher σ, as compared to the homogeneous SSS case.

Robustness.
In the previous Sections we have addressed the dependence of relay synchronization in a multiplex on the dynamical and structural layer heterogeneity, and proved that the phenomenon still holds even when the intermediate layer has a completely different structure and dynamics than the mirrored ones. The present section is devoted instead to assess the robustness of relay synchronization against structural changes by means of a de-multiplexing process of the layers, that is, against performing a progressively shutting down of the inter-layer links such that a fraction of nodes in each layer is not linked to their counterparts in the other layers. In addition, we also investigated several cases to test the robustness against mismatches in the oscillators parameters to closely resemble real experimental conditions. Structural robustness. To closely check this process, we initially consider a 3-layer multiplex with identical topology (EEE or SSS). We choose the inter-and intra-layer couplings to guarantee a relay synchronous state with the layers fully multiplexed. Then, we proceed to disconnect one by one the inter-layer links according to the nodes degree ranking, both in the ascending and the descending order, and re-evaluate in every step the state of the relay synchronization by measuring the E −1,1 error. Figure 6(a) reports the averaged evolution of E −1,1 as a function of the number of multiplexed nodes, after having performed the whole de-multiplexing process for 10 different network realizations. It can be seen that, starting from a situation with E −1,1 = 0, the EEE multiplex configuration (blue void markers) loses the synchronization immediately with just a few of inter-layer links being removed. On  the other hand, relay synchronization is resilient in SSS triplex configurations also when more than 30% of the nodes are not multiplexed.
A more detailed view can be obtained from Fig. 6(b), where the number of multiplexed nodes needed to support the relay synchronization is represented as a function of the intra-layer coupling σ 0 = σ 1 . As expected, when the coupling is weak, all the N nodes need to be linked to preserve relay synchronization. However, as the interaction within the layers increases, the intra-layer connectivity helps to maintain a synchronous state despite an increasing number of nodes are being de-multiplexed without damaging the coherence between the outer layers. In Fig. 6(b), we can see that for both the EEE (blue void markers) and the SSS (black full markers) triplex configurations, removing the links between layers connecting nodes with higher degree (descending degree ranking, circle markers) is much more robust than following an ascending degree ranking (square markers). This is indeed a very interesting result: relay synchronization in a multiplex network is supported by the low degree nodes, while the hubs can be safely disconnected without perturbing the transmission. This is notably evidenced in the SSS case (black full squares) where after having removed the 40% of the inter-layer links connecting the highest degree nodes, the relay synchronization is still supported by the multiplex structure connected through the lowest degree nodes.
Once we have singled out that, from our trial rankings, the descending degree ranking is the most convenient way to de-multiplex part of the network without loosing coherence, we proceed our study by evaluating the impact of having a relay layer with different topology from the outer layers, as we did in the previous Section 4. In this scenario, we have two possible descending degree rankings, the one dictated by the relay layer and the one dictated by the outer layers. The results are summarized in Fig. 7 where we plot, as in Fig. 6(b), the number of nodes that need to be linked to maintain synchronization as a function of σ 0 = σ 1 . For the sake of comparison, we added the curves for the homogeneous EEE and SSS (full markers) multiplex configurations, together with the data for the mixed ESE and SES (void markers) layers. Notice that the chosen inter-layer coupling λ = 0.23 is well above threshold for all the cases, as it can be derived from Fig. 5. All the reported evidence indicates that the introduction of a relay layer with a topology different from that of the outer layers has little influence on the minimum number needed to support the relay synchronization, as long as the first removed inter-layer connections correspond to the hubs in the outer layers (blue and back curves). Curiously, the alternative of using the relay layer topology to rank the degree of the nodes, destroys the coherence between the outer layers as soon as a tiny fraction of links is removed (red curves). Therefore, the relay synchronization in a multiplex is very unstable if just a few links connecting nodes which are hubs in the relay layer are removed. Notice that this unlinking criterion is equivalent to randomly disconnect the multiplex. Therefore, the robustness of the relay synchrony relies mainly in the low degree nodes of the external layers. The relevance of the low degree nodes in controlling the dynamics of complex networks has been pointed out in other contexts 43,44 .
Dynamical robustness. In order to explore a more realistic situation of a multiplex composed of non identical oscillators, we have taken further the generalization by introducing some heterogeneity in the node's parameter a k , that is, a k = a k,i . The node values of the outer layers's were randomly picked from the interval a ±1,i = [0.20, 0.21].
We prepared three different setups depending on how the parameter setting is introduced in the outer layers: i) a 1,i = a −1,i but ≠ a a i j 1, 1, , that is, oscillators are not identical within each external layer but the two external layers are symmetric; ii) a 1,i = a 1,j but a 1,i ≠ a −1,i , that is, layers are composed of identical oscillators but there is a mismatch between oscillators in layer k = 1 and their corresponding replicas in layer k = −1; and iii) SCIeNtIfIC REpoRtS | (2018) 8:8629 | DOI:10.1038/s41598-018-26945-w a 1,i ≠ a 1,j ≠ a −1,i , that is, combining cases i) and ii), a fully random case, not preserving the symmetry neither the parameter homogeneity in each layer.
Results are summarized in Fig. 8, where the synchronization error between the outer layers is plotted vs the inter-layer coupling λ. For comparison, the perfectly homogenous case with all layers having identical oscillators (red circles) is also included. In the inset, a zoom of the main figure, it is clear that breaking the symmetry between the outer layers slightly deteriorates the synchronization, as it occurs in cases i) (green squares) and iii) (purple diamonds). However, introducing some heterogeneity in the a k parameter (full cyan circles) in the layer's oscillators does not modify the synchronization threshold at all, as long as each oscillator in one outer layer is identical to its replica in the mirrored layer. We will actually show that our experimental results can be framed within this latter case.

Experimental validation
Finally, we report experimental evidence of relay synchronization in a multiplex of nonlinear electronic circuits, with the setup sketched in Fig. 9 (left). The array is made of 21 Rössler-like circuits arranged in three layers of 7 nodes, with the relay layer having different topology as the outer layers. Each layer has two different electronic couplers, one for the coupling among nodes in the same layer (σ e ) and the second for the interaction of each node with its replica in the other layers (λ e ). The chaotic dynamics of the circuits is well approximated by the three variables (v 1 , v 2 , v 3 ) obeying 23 :   where G v i 1 is a nonlinear gain function given by: where the parameter values are gathered in Table 1. The reader is referred to ref. 45    the PC for further analysis. Next, the coupling between the inter-layer (λ e ) increases one step (0.01), digital pulses are sent to the potentiometers corresponding to that coupling and decreases the resistance in 100Ω each time it passes through this state, until the maximum value of λ e is reached (Ω in potentiometers). When the entire run is finished, σ e is increased by one step, and another cycle of λ is initiated. The whole procedure is repeated until each coupling reached its maximum value. The experiment is controlled from a PC with the LabVIEW software.
The experimental results are summarized in Fig. 10. The top panels represent the averaged experimental inter-layer synchronization error for the outer layers E −1,1 (left) and between the relay and one of the outer layers E 0,1 (right), for all the experimental range of intra-layer σ e = [0, 0.6] and inter-layers λ e = [0, 0.6] couplings. Even though the system is unavoidably affected by noise and parameter mismatch in the electronic components, for high enough λ e the value of E −1,1 is well below E 0,1 and therefore the inter-layer relay synchronization is verified in our experimental setup. Superimposed to the colormaps, we also have drawn the isoline for E = 0.12 in both panels (white lines), showing that the threshold λ ⁎ e value for which E −1,1 and E 0,1 are below the value of the isoline is always smaller in the E −1,1 case. The fact that the perfect synchronization between the two outer layers is never achieved agrees with our numerical predictions reported in the dynamical robustness section and cleary visualized in Fig. 8.
For a clearer view, in the bottom left panel we have just plotted E −1,1 and E 0,1 as a function of λ e for a fixed intra-layer coupling σ e = 0.5 (corresponding to the blue and black dashed lines in the respective colormap panels in the upper part of Fig. 10), showing that E −1,1 monotonically goes to zero and is always below E 0,1 .
Finally, in the bottom-right panel of Fig. 10 we plot both errors, E −1,1 and E 0,1 , as a function of σ e for a fixed value of the intra-layer coupling λ e = 0.5 (vertical cuts in red and magenta in the colormap plots). That is done in order to show the effect of increasing the interaction in the intra-layer connectivity. Similarly to what observed in Fig. 5, promoting the topological difference between layers as σ e increases rises the synchronization threshold.

Discussion
The synchronous behavior of groups of units in a complex system is often a signature of a common functional involvement. Theoretically, this has been studied in the framework of modular networks 47 , where nodes densely connected among them in the same mesoscale structure usually share similar functions [48][49][50][51][52][53][54][55] , or in the context of cluster synchronization associated to the existence of network symmetries and where modularity is not relevant 35,56 . However, most of these studies disregard the possible different nature of the involved nodes and links, a common feature in real systems. For instance, in the brain coherence is observed involving different cell types and electrical/chemical synapses.
In addition, cluster and modular synchronization requires the direct connection between the synchronizing nodes. However, long distance coherence between complex mirrored structures mediated through non-synchronous differentiated ones plays a key role in the functioning of several real-world systems. Zero-lag synchronization has been indeed observed between distant areas of the cortex 37,38 , with the thalamus acting as a coordinator, and the transcendental role of symmetry in its dynamics has been lately pointed out in other contexts like in the evolution of complex developmental systems 57,58 .
In this work, we have overcame these limitations by using a multilayer description in the study of distant synchronization in heterogeneous ensembles. Within this framework we have accounted for, besides considering different coupling functions between the dynamical units, the impact of having topologically different layers and heterogeneity in the node dynamics. We have implemented the concept of relay synchronization to the case of a multiplex network, showing that the intermediation of a relay layer can lead to inter-layer synchronization of a set of paired layers, both topologically and dynamically different from the transmitter. The phenomenon can be extended to indefinitely higher order relay configurations, provided a mirror symmetry is preserved in the multiplex. The coherent state is very robust to changes in the dynamics, topology, and even to strong multiplex disconnection. In this latter scenario, we show that the lower degree nodes in the synchronized outer layers are responsible for resilience of the synchronous state, while hubs can be safely de-mutiplexed. Finally, we experimentally validated our results in a multiplex network of three layers of electronic oscillators. Our results provide a new path for starting the study of the role of symmetries in setting long distance coherence in real systems.