Relay Synchronization in a Weighted Triplex Network

: Relay synchronization in multi-layer networks implies inter-layer synchronization between two indirectly connected layers through a relay layer. In this work, we study the relay synchronization in a three-layer multiplex network by introducing degree-based weighting mechanisms. The mechanism of within-layer connectivity may be hubs-repelling or hubs-attracting whenever low-degree or high-degree nodes receive strong inﬂuence. We adjust the remote layers to hubs-attracting coupling, whereas the relay layer may be unweighted, hubs-repelling, or hubs-attracting network. We establish that relay synchronization is improved when the relay layer is hubs-repelling compared to the other cases. We determine analytically necessary stability conditions of relay synchronization state using the master stability function approach. Finally, we explore the relation between synchronization and the topological property of the relay layer. We ﬁnd that a higher clustering coefﬁcient hinders synchronizability, and vice versa. We also look into the intra-layer synchronization in the proposed weighted triplex network and establish that intra-layer synchronization occurs in a wider range when relay layer is hubs-attracting.


Introduction
Over the last few decades, the study of complex networks has been a topic of great interest due to its application in diverse scientific fields [1,2]. Recent advances in network science revealed that many systems in nature could be represented by multi-layer networks [3,4], where elements are interconnected through various links between more than one layer. Such a framework is well known for illustrating many real-world systems, such as social networks [5], transportation networks [6], power-grid networks [7], ecological networks [8], and neuronal networks of the brain [9,10]. A particular subclass of multi-layer networks are multiplex networks where each layer contains the same number of nodes, and one-to-one inter-layer connections between the replica nodes from neighboring layers are allowed.
Relay synchronization in a multiplex network has received much attention from researchers in the past few years. This synchronization phenomenon occurs between two distant layers, not directly connected but only through a relay layer. First reported in a star-like network [29,30], this type of synchronization has also been detected in laser systems [31], circuits of chaotic oscillators [32], and between remote pairs of nodes within the network [33]. Besides, it is of great importance in the human brain networks, where the thalamus [34] acts as a relay between distant cortical areas. Recently, relay synchronization has been illustrated in multi-layer and multiplex frameworks with different intra-layer topologies, such as Erdős-Rényi and scale-free [35]. In Reference [36], random inhomogeneities of the small-world type have been introduced in the network layers to scrutinize the robustness of relay synchronization. Continuous and discrete-time systems have been used as individual units to induce spatiotemporal structures in triplex networks [37,38]. However, in the above cases, the network topology of both the relay and remote layers are represented by unweighted networks. Although weighted networks, where coupling weights are heterogeneous, are present in many real-world networks, such as the network of scientific collaborations [39], structural and functional networks of the brain [40], transportation networks [41], epidemic spreading [42], ecological networks [43], etc., we would like to emphasize that the recent progress in this direction has been achieved considering weighted mono-layer networks [44][45][46][47][48] and weighted duplex networks [49][50][51]. However, the effect of weighted networks upon synchronization in the multiplex network with more than two layers is less studied. So, the lack of works on the synchronization in weighted multi-layer networks encouraged us to unveil the effect of weighted networks upon the relay synchronization in a multiplex network.
In this regard, we would like to highlight the very recent works by Estrada et al. [52,53], in which the authors have introduced two opposite approaches to explore the degree-biased mechanisms on the network dynamics based on the generalization of the graph Laplacian operator. The Laplacians are defined so that the pairwise weights of connectivity matrix are introduced as a ratio of degrees of two coupled nodes. As a result, an unweighted graph is transformed into a weighted directed one following two distinct methods. One forces the network to obey high-degree nodes, and the other is biased towards the low-degree nodes. In the first case, referred as hubs-attracting Laplacian [52], a high-degree unit receives weaker influence from its low-degree neighbor than from a high-degree one. Opposed to the above, the low-degree units produce more potent effects in the second case, illustrated by the hubs-repelling Laplacian [53]. These findings have motivated us to dig into the comparison between the weighted directed networks generated by different degree-based weighting mechanisms to achieve synchronization in a relay multiplex network.
In this letter, we consider a three-layer multiplex (triplex) network with intra-layer coupling set up by a random Erdős-Rényi (ER) topology with individual nodes as chaotic Rössler oscillator. The intra-layer connectivity is set up in such a way that the corresponding Laplacian matrices are equal to either hubs-attracting or the hubs-repelling Laplacian matrices. In particular, we set the remote layers to be hubs-attracting and consider three configurations of relay layer: hubs-attracting, hubs-repelling, and unweighted network. Further, we estimate the onset of relay synchronization for all three configurations. Our study shows that the relay synchronization is improved when the relay layer is hubsrepelling compared to the other two cases. Using analytical tools, we demonstrate the necessary conditions for the synchronized state. We also explain the enhancing behavior in terms of clustering coefficient analysis of the relay layer. Apart from this, we also examine the emergence of intra-layer synchrony for all three configurations.

Mathematical Model
We consider a triplex network, schematically presented in Figure 1. The layers are denoted by k = −1, 0, 1, where k = ±1 indicates remote layers and k = 0 corresponds to the relay layer. Each layer contains N = 100 nodes following m-dimensional identical dynamical systems. The nodes are mutually connected with their corresponding replica nodes from the neighboring layers. The states of the kth layer are represented by the vectors X k = (x k,1 , x k,2 , . . . , x k,N ), where x k,i ∈ R m and i = 1, 2, . . . , N. The dynamics of the ith node in each layer is governed by the system of ordinary differential equations: represent the autonomous evolution of the uncoupled oscillator, intra-layer, and inter-layer coupling matrices, respectively. In the current study, we suggest that the autonomous evolution of each node in the triplex network (1) obeys the dynamics of chaotic Rössler oscillators, for which F(x k,i ) is given by (2) In system (2), the control parameters are fixed and set equal for all oscillators as a = b = 0.2, c = 5.7, which provide chaotic dynamics. The intra-layer coupling matrix G for all the layers is taken as linear diffusive coupling through x-variable, i.e, G = [1, 0, 0] tr , where tr denotes the transpose of a matrix. The inter-layer coupling is considered to be through all the state variables. Thus, the inter-layer coupling matrix H is an identity matrix. Strengths of intra-layer and inter-layer coupling are introduced as σ and λ, respectively.
Laplacian matrix of the k th layer is defined as L ij . The parameter β determines if the layer remains unweighted network (U, specifies that the high-degree nodes produce strong influence on their low-degree neighbors, and vice versa.
ij defines the opposite mechanism, the low-degree units have a strong impact on their high-degree neighbors, and vice versa. The Laplacian obtained for β = ±1, respectively, corresponds to hubs-repelling [53] and hubs-attracting [52] Laplacian matrices, and associated networks are called hubs-repelling and hubs-attracting networks.
Here, the connection topology of each layer is given by Erdős-Rényi (ER) random graph with the edge connectivity probability p rand = 0.05. We transform the network of each layer to a directed weighted graph by introducing the mechanism of hubs-attracting (A) and hubs-repelling (R) Laplacian. The remote layers are transformed to the hubsattracting network, and the middle (relay) layer could be transformed to one of three configurations: Hubs-repelling (R, β = −1), unweighted (U, β = 0), and Hubs-attracting (A, β = 1). We classify the layer stacks regarding the weighting sequence of each layer. For instance, the triplex network with remote layers Hubs-attracting and middle layer Hubs-repelling (β = −1) will be denoted as ARA. Consequently, the triplex networks with middle layers unweighted and Hubs-attracting will be denoted as AUA and AAA, respectively. We will scrutinize the relay synchronization for all three configurations: AAA, AUA, and ARA in the next section.

Valuation of Relay Synchronization Error
To analyze the relay synchronization in our networks, we employ the synchronization error E 1,−1 as: where · denotes the Euclidean norm, t tr determines the transient of the numerical simulation, and T is a sufficiently large positive number. To measure the synchronization errors, the time interval is taken over 1 × 10 5 units after an initial transient of 4 × 10 5 units, i.e., t tr = 4000 and T = 1000. In all our simulations, we consider the threshold value of synchronization error to be 10 −5 for synchrony. The triplex network (1) is solved numerically using the Runge-Kutta-Fehlberg algorithm with integration step dt = 0.01 and fixed random initial conditions chosen from the interval [−1, 1]. For each result, 20 network realizations with random initial conditions are considered.
We first scrutinize synchronization error for a fixed value of intra-layer coupling strength σ = 0.02 by varying inter-layer coupling λ. Figure 2 displays three curves for three different configurations: AAA, AUA, and ARA, respectively. For the AAA configuration, red dotted line, the synchrony is achieved at the critical coupling λ ≈ 0.074. For the AUA configuration, black dotted curve, the critical value shifts towards a slightly lower value of λ. At the same time, for the ARA case, blue dotted curve, the critical value of synchronization decreases further -the synchrony is reached at λ ≈ 0.07. Thus, one can conclude that, for a fixed value of intra-layer coupling, the relay synchronization enhances as the value of β decreases from 1 to −1 within the relay layer, and the remote layers are fixed to be hubs-attracting (A) network.
Therefore, without introducing parameter mismatch in the relay layer or introducing a time delay between the layers, changing the weighting mechanism between the middle (relay) layer nodes may effectively improve the relay synchronization in multiplex networks.

Stability Analysis of Relay Synchronization
Thereafter, we analytically derive the conditions for stability of the relay synchronization using Master Stability Function [54] approach. The synchronization solution is given by, X 1 = X 0 .
Let δX(t) = X 1 (t) − X −1 (t) be the vector describing the difference between the dynamics of the mirror layers. Considering δX = {δx 1 , δx 2 , . . . , δx N } to be a small perturbation of synchronous solution, we have linearized Equation (1) as: where J is the Jacobian operator, and X = {x i } represents the synchronous state In the above equation, x 0,i is the state variable of the i th node in the relay layer, which satisfies the equation of motioṅ In the variational Equation (5), the state variable evolves transverse to the synchronization manifold. Therefore, the Lyapunov exponents of Equation (5) depict the stability of the synchronization manifold. The maximum Lyapunov exponent (MLE) as a function of the parameters σ, λ gives a necessary condition for the stability of the inter-layer synchronization solution. For a stable solution, perturbation along all the transverse directions must die out, i.e., the values of MLE should be negative. Solving the linearized Equation (5) along with the non-linear Equations (6) and (7) gives the MLEs transverse to the synchronization manifold. In Figure 3, we have plotted the MLEs for the same parameter values as in Figure 1 to validate the stability condition. It clearly shows that the MLE passes zero-level at the exact value of λ where E 1,−1 becomes less than 10 −5 for each of the three configurations considered. Therefore, the analytical stability condition suitably agreed with our numerical results.

Delineation of Relay Synchronization in Parameter Space
To better understand the enhancement in the relay synchronization, we have plotted the variation of synchronization error E 1,−1 in the (λ, σ)-plane for all the three network configurations. The results depicted in Figure 4 show that the region of synchronization increases in ARA (Figure 4c) configuration compared to the other configurations, AAA ( Figure 4a) and AUA (Figure 4b). Thus, rigorously plotting the synchronization error, we can perceive the synchronization scenario for all three configurations. It reaffirms our statement about an enhancement of relay synchronization in ARA configuration of the multiplex network as compared to AAA and AUA configurations.

Structural Analysis
Apart from quantifying the synchronization error and linear stability analysis, for a profound understanding of enhancement in the relay synchronization, we also go through the structural analysis of the middle (relay) layer for all three cases. With increasing β from −1 to 1, we analyze the average clustering coefficient, CC = 1/N ∑ N i=1 C i , where C i is a local clustering coefficient, within the relay layer. As for non-zero β, the network becomes a weighted directed network, and we follow the formula for local clustering coefficient for weighted directed network proposed in Reference [55]: where d in i and s in i are the in-degree and input strength of ith node, respectively. Here, we have taken only the in-component of clustering coefficient, as the hubs-attracting network's input strength corresponds to the hubs-repelling network's output strength, and vice versa. In Figure 5, we plotted the average clustering coefficient for different values of β within the relay layer. The largest value (CC = 0.0649) of clustering coefficient corresponds to β = 1 when the relay layer is Hubs-attracting, and for β = −1, when relay layer is Hubs repelling, the value of clustering coefficient is lowest (CC = 0.0630). It is a well-known fact that the clustering coefficient plays a crucial role in synchronization. Larger values of clustering coefficient hinder the synchronization, and smaller values determine the improvement [56]. Due to large clustering, the network splits into dynamical clusters that follow different trajectories and, as a result, impedes synchronization. Hence, the result depicted in Figure 5 can reasonably explain improving the relay synchronization with ARA configuration.

Emergence of Intra-Layer Synchronization State
Now, we investigate the intra-layer synchronization in the multiplex network (1). For intra-layer synchronization, each individual layer evolves on the same time evolution, which occurs when each individual oscillator in each layer of the multiplex network is appropriately coupled [17,57]. To quantify the intra-layer synchronization state, we introduce the synchronization error as, is the synchronization error of each layer k. E intra ≈ 0 indicates the emergence of intra-layer synchronization in the triplex network.
We have plotted the intra-layer synchronization error E intra in Figure 6 by varying the intra-layer coupling strength σ for the three configurations: AAA (red dotted line), AUA (black dotted line) and ARA (blue dotted line), respectively. From the figure, we can elucidate that, for all the three cases, the triplex network shows two transitions-one from desynchrony to synchrony and other one from synchrony to desynchrony state. So, the intra-layer synchronization appears in a bounded range of intra-layer coupling strength σ. For AAA configuration, the first transition occurs at intra-layer coupling σ = 0.13 and the second transition takes place at σ = 0.44. For AUA configuration, the intra-layer synchronization occurs in the range [0.17, 0.38] of intra-layer coupling strength. The range shrinks further for the ARA configuration for which synchronization occurs in [0.19, 0.26]. One point of note is that the intra-layer synchronization occurred simultaneously [58], and desynchronization to synchronization also occurred at the same value of intra-layer coupling strength σ.
Therefore, intra-layer synchronization for the multiplex network emerges in a bounded range of coupling strength. The range is wider for AAA configuration compared to AUA and ARA configurations.

Discussion and Conclusions
In summary, we have explicitly addressed the emergence of relay synchronization in weighted triplex networks of chaotic non-linear oscillators. We constructed withinlayer connection topologies following the Erdős-Rényi algorithm for N = 100 nodes. The pairwise intra-layer interactions are scaled such that the influence produced by the ith node on the jth node is proportional to the ratio of their degrees The current study evidences an enhancement in the relay synchronization if the middle layer exhibits hubs-repelling interactions compared to hubs-attracting and original unweighted couplings. We have supported numerical findings by analytical treatment of the synchronized state stability conducted via master stability function. Moreover, we have also scrutinized the topological properties of the intra-layer connection within the relay layers, which helps in understanding that, in the case of hubs-repelling interactions, the middle layer exhibits a lower clustering coefficient. The latter contributes to a more efficient relay synchronization in a considered triplex network. Lastly, we investigated the intra-layer synchronization for the multiplex network which emerges in a bounded range of intra-layer coupling strength and the range of synchrony expands if the relay layer is hubs-attracting compared to unweighted and hubs-repelling couplings.
We expect that our investigation considering such mechanisms of network weighting as hubs-attracting and hubs-repelling will open the door for a new perception of various collective behaviors emerging in a multi-layer network. Data Availability Statement: All data needed to evaluate the conclusions in the paper are present in the paper itself. Additional data related to this paper may be requested to the corresponding authors.