Interlayer Hebbian Plasticity Induces First-Order Transition in Multiplex Networks

Adaptation plays a pivotal role in the evolution of natural and artificial complex systems, and in the determination of their functionality. Here, we investigate the impact of adaptive inter-layer processes on intra-layer synchronization in multiplex networks. The considered adaptation mechanism is governed by a Hebbian learning rule, i.e., the link weight between a pair of interconnected nodes is enhanced if the two nodes are in phase. Such adaptive coupling induces an irreversible first-order transition route to synchronization accompanied with a hysteresis. We provide rigorous analytic predictions of the critical coupling strengths for the onset of synchronization and de-synchronization, and verify all our theoretical predictions by means of extensive numerical simulations.


Introduction
The study of synchronization, a collective motion of initially un-identical units interacting on network structures has provided a deeper understanding on the underlying processes (and the nature of transition) taking place on a wide range of physical and biological systems [1]. Second-order-like transitions in systems of phase oscillators are arXiv:2010.09424v1 [nlin.AO] 19 Oct 2020 frequent in nature, whereas abrupt, discontinuous, first-order-like transitions are not so common. However, a plethora of recent studies has established that first-order-like transitions can be achieved in networked oscillators under certain circumstances inducing frustration mechanisms preventing synchronization between connected oscillators, and thus blocking the formation of giant clusters during the transition. Such a transition (called explosive synchronization) usually involves an abrupt formation of the giant synchronized cluster, and then its irreversible abrupt de-synchronization, yielding a hysteresis loop. The hypersensitivity or abruptness arisen due to a minuscule change in the interaction strength makes this process unmanageable and calamitous in many circumstances. Examples are, for instance, blackouts in the power grid [2], breakdown of the internet [3], episodes of Fibromyalgia chronic pain [4] and epileptic seizures [5] in the human brain. The explosive synchronization (ES) transition is shown to emerge in networked oscillators with microscopic correlation between their frequency and network's structural property such as degree [6] and coupling strength [7], or by inclusion of inertia [8,9,10] and noise [10,11], or by the presence of a fraction of adaptively coupled oscillators [12,13,14,15,16,17], mean-field diffusion [18], traffic processes [19], the presence of nearest-neighbor competitive interaction or symmetrybreaking interaction [20] and anti-Hebbian adaptation of link weight [21].
Despite being a very useful framework to investigate phenomena such as synchronization and percolation, isolated networks are often unable to precisely represent the behavior of complex systems involving different types of interactions among the same set of interacting elements. Multilayer or multiplex networks are the right candidates to map such systems, as they comprise different types of connections or processes among a common set of nodes, interconnected through different layers [22,23,24,25,26]. The advent of multiplex network has made it possible to investigate the impact of one type of process (layer) on other interdependent processes (other layers) such as synchronization, percolation, epidemic spreading, etc. In the same breath, a few techniques have been employed inducing ES in one or all dynamical layers. For instance, ES is common in adaptively coupled interdependent layers [12], in random walker dynamics [27], in presence of a delay [28] or interlayer adaptation through order parameter of the layers [29].
The mechanism of adaptation plays a crucial role in the development and function of many natural and artificial systems. In neuroscience, most of the experimental findings suggest that synaptic adaptation between neurons is the basis of learning and long-term memory [30,31]. First proposed by Hebb [32] and later supported by experimental evidences [31,33,34], such adaptation consists in the fact that the synaptic coupling between two neurons is strengthened or weakened if the two neurons are simultaneously firing or not depending on the pinpoint relative timing of presynaptic and postsynaptic spikes. If the relative spike timing is coded in terms of the phases of the oscillators, a neural network then can be delineated as a network of phase oscillators. Thus, the adaptive evolution of the neural network occurs through the alterations of synaptic connections between neurons. The adaptation of connection-weights following such Hebbian learning mechanism leads to the occurrence of cluster states [35,36,37] or meso-scale structures [38,39] underlying synchronization in complex networks.
In this paper, we investigate a multiplex framework inspired by such Hebbian adaptive learning rule. In our network, the weights of the links between interconnected layers are adaptive and governed by the instantaneous phase-difference of the interconnected nodes. Strikingly, the Hebbian learning rule divides the population of interlayer links into weights and anti-weights for low intra-layer interaction strength. The existence of the inhibitory interlayer weights drives the multiplexed layers adopt the first-order transition route to synchronization accompanied with a hysteresis. The Hebbian learning mechanism also provides a great amount of control over the abruptness and the width of associated hysteresis by means of learning parameters. It is further revealed that the critical coupling strength for the onset of synchronization does show explicit dependence on the learning parameter while that for the onset of desynchronization remains independent of it. The numerical assessments of both the critical coupling strength have shown a good match with their respective analytical predictions. The proposed recipe based on Hebbian adaptation is capable of triggering first-order transitions in all homogeneous networks.

Model
Let us start with investigating how the Hebb's neural learning mechanism governing the inter-layer weight (strength) affects the phase transition in the multiplexed layers. In order to achieve this, we consider a multiplex network composed of two layers of the same size N . The dynamics of nodes in each layer is governed by the Kuramoto model [40]. The inter-layer link weight D [ii] between node i in a layer and its counterpart in the other layer is adaptive in nature. Hence, the evolution of the phase oscillators is governed byθ where i = 1, ..., N , subscripts 1 and 2 stand for the two distinct layers, and λ represents intra-layer coupling strength, here λ 1 = λ 2 = λ. The instantaneous phase of the i th node is denoted by θ l follows uniform or unimodal distribution g(ω l ). The intra-layer connectivity between the nodes following a network topology is encoded in the adjacency matrix A l such that A [ij] l = 1 (0) if i th and j th nodes are connected (disconnected). The dynamically adaptive weight D [ii] of an inter-layer link between each pair of interconnected (mirror) nodes in the two layers is determined by the following Hebbian learning rulė where α ∈ [0, 1] is a factor which amplifies the amount of learning if the two nodes  Figure 1. Schematic representation of the duplex network (l = 1, 2) having non-adaptive or adaptive inter-layer coupling between interconnected nodes. The nonadaptive case (learning rate ε = 0) leads to a continuous transition, whereas the Hebbian learning rule (ε = 0) yields a discontinuous route to synchronization and then an irreversible route to desynchronization. r f and r b denote order parameter in the forward and backward continuation of coupling strength λ, respectively.
are synchronized and ε ∈ [0, 1] is the learning rate. The D [ii] in Eq.2 prevents the inter-layer coupling weight from increasing or decreasing without bounds. Hence, the supra-adjacency matrix of the multiplex network is adaptive, and includes un-weighted intra-layer links and adaptive weighted inter-layer links: where I is the identity matrix.
To track the level of coherence in the system, we define the global order parameter r l for a layer l in terms of the average phase ψ l as (4) r = 1 represents a completely synchronous state, while r = 0 implies total incoherence. In a similar way, one can define in-phase (one-cluster) global order parameter for the entire multiplex network as Furthermore, to capture the degree of anti-phase (two-cluster) global synchronization, a new global order parameter R 2 (dipole moment of the distribution of anti-phases) [35,41] is defined as follows 0.00 0.05 0.10 0.0 0.5 where Here, R measures both the one-cluster and the two-cluster synchronization. Hence in order to determine the degree of two-cluster synchronization, the term R is adjusted in Eq.6 for one-cluster synchronization.

Results
We consider a multiplex network made of two Erdös-Rényi (ER) random networks [42] having average intra-layer connectivity k 1 = k 2 = 10. The nodes in both layers are assigned initial phases and natural frequencies drawn from a uniform random distribution such that is the number of inter-layer links a node in one layer can have with the nodes in other layer. Here we adopt the simplest form of a multilayer network, i.e. d [i] = 1, ∀i.
First, we investigate the impact of the learning rate ε on intra-layer synchronization in the multiplexed layers. Fig.2 illustrates the behavior of the order parameter (for both forward and backward transitions) for the two layers as a function of coupling strength λ for different values of learning rate ε sustained with a fixed choice of learning factor α. It is found that the two layers undergo a first-order transition (ES) for different values of ε. It is also unveiled that the backward critical coupling strength λ b c is independent of the rate ε. However, the forward critical coupling strength λ f c at first decreases with the increase in ε, then no further change is observed for higher values (0.5 − 1) of α. Hence slower learning rates ε yield slightly wider hystereses.
Next, the impact of the amplification factor α with a fixed choice of ε is reported in Fig.3. The two layers follow a first-order transition with significantly different hysteresis width for different values of α. At very low α, a second-order transition is observed. With the increase in α, the forward critical λ f c significantly increases while the backward critical coupling λ b c remains independent on α as well. Thus, each increase in α leads to a broader hysteresis width.
To gather information at a microscopic level on the origin of the first-order transition, we investigate the behavior of 1 |, in phase (R 1 ) and anti-phase (R 2 ) global order parameter in Fig.4. The first column of Fig.4 reports the final state of D [ii] in the forward continuation of λ for different values of α at a fixed rate ε. It is apparent that the stationary values of In the incoherent state (λ < λ f,b c ), the inter-layer link population is mainly divided into two clusters of anti-weights, namely −α and α. Besides, there exists a neutral cluster around D [ii] → 0 containing a small fraction of link population, which increases with the decrease in the value of α. Interestingly, D [ii] and −D [ii] populations are almost equal in size. On the contrary, the inter-layer link population tends to converge to −α in the coherent state (λ > λ f,b c ). The inter-layer link population having D [ii] → −α in the incoherent state originates a frustration at the respective interconnected end nodes. The triggered frustration at the end nodes of the inhibited (subjected to −α) inter-layer links curbs the formation of the largest synchronous cluster in their respective layers until a threshold for λ is reached. The larger the magnitude −α, the stronger the triggered frustration. Since a large value of −α imparts a stronger inhibition during the forward continuation of λ, a stronger and stronger λ is required for the abrupt formation of the largest synchronous cluster with each increase in α. Thus, the onset of first-order transition is witnessed at a larger forward critical coupling λ f c with each increase in α. In the second column of Fig.4, we study the distribution P (∆θ 1 |, the difference between phases of the interconnected nodes for different values of α while keeping ε fixed. It unveils that in the incoherent state (λ < λ f c ) belonging to an intermediate or higher value of α, two phase-clusters in each layer exist: one corresponding to θ 1 and the other to |θ In addition, the neutral cluster population remains distributed between these two clusters, whose population increases with the decrease in α. Hence, each multiplexed layer stays in a bi-clusters state in the incoherent state. Nevertheless, the P (θ [i] ) in the coherent state (λ > λ f c ) unveils a sharp single peaked (unimodal) distribution centered at |θ 1 | = π for any value of α. It implies that both layers are locked to two different phases at a mutual difference of π, i.e., the multiplex network comprises two anti-phase layers in the coherent state.
The third column of Fig.4 illustrates R 1 (Eq.5) and R 2 (Eq.6) for different values of α. In the incoherent state one has that R 1 → 0 and R 2 → 0, as there exists two antiphase clusters with θ [i] = 0 and θ [i] = π in each layer. In the coherent state, however, there exists a single cluster centered at θ [i] = π, hence still one ha R 1 → 0, while R 2 displays an abrupt jump at the critical coupling strength. The height of the abrupt jump for R 2 increases as α increases. Now, sinceḊ [ii] = 0 (from Eq.2) in the stationary state, one obtains (for ε = 0) Further, it is revealed from Fig. 4  1 |= arccos(±1) = 0, π. Hence, the existence of the inhibitory D [ii] −α gives rise to anti-phase mirror-populations in the two layers in both incoherent and coherent states.
Phase Diagrams: The data shown in Fig.2 highlight fact that any ε = 0 is capable of inducing a first-order transition in the system and there exists a marginal difference in the critical coupling strength λ f c for different values of ε. Fig.5 reports on how the hysteresis width |r f l − r b l | varies in the α − λ space for different values of ε and for multiplexes composed of two homogeneous (ER-ER) and by one homogeneous (ER) and one heterogeneous (Barabási-Albert (BA) [43]) networks. The ER-ER configuration exhibits a first-order transition for intermediate and higher values of α. Also, the associated hysteresis width for both the layers increases with an increase in α. Besides, a slower rate ε yields a wider hysteresis width for a given value of α as compared to the faster rate ε. For the ER-BA configuration, a rather weak hysteresis is observed for the ER layer in the range α = 0.5 − 1, while the BA layer does not feature a first-order transition route to synchronization.

Analytical treatment
To obtain an analytical expression for the order parameter, we take into account, for simplicity, a multiplex network consisting of two globally-connected (GC) layers l ∈ [1,2] so that A ij l = 1/N in model Eq.1. The intrinsic frequencies of the nodes in both GC layers are selected from a uniform distribution The order parameter for each layer l ∈ [1, 2] is then defined as Now Eq.(1) can be rewritten in the mean-field form using Eq.(8) aṡ In the stationary stateḊ [ii] = 0, for ε = 0, thereby Eq.2 leads to Hence, evolution of the nodes in the stationary state is ruled by the following selfconsistent equationṡ Next, we gather from numerical simulations (see Fig.5) that the distribution 1 | for the two GC layers in the coherent state (λ > λ b c ) follows a peaked (unimodal) distribution with its mean at 0 and standard deviation ∆θ(λ). Therefore, the inter-layer term in the stationary state can be expressed as σ = α 2 sin(2∆θ), and the model Eqs.9 can be rewritten aṡ θ In this way, stationary σ = σ(α, λ) accounts for the maximum possible inter-layer contribution in the evolution of phases in either layers. When the phases in each layer are locked to their respective mean-fieldsθ , which leads to the following set of the conditions (referred as CondI) to be satisfied simultaneously by the locked oscillators contributing to r l ; CondI: The order parameter in Eq.8 for the phase-locked oscillators can be expressed as which can be further simplified as Now, ∆θ → 0 (∆θ 0.002 − 0.003) for any λ > λ b lc corresponding to locked state, hence σ → 0 (see Fig.6).
Backward Critical Coupling: In the continuum limit N → ∞, the order parameter can be rewritten in its integral form Interlayer Hebbian Plasticity Induces First-Order Transition in Multiplex Networks 11 For a uniform distribution g(ω l ) = 1 2γ for |ω l | < γ and since Ω l → 0, the order parameter in Eq.15 reduces to If ∆θ [i] ∼ π for the pairs of oscillators locked in their respective layers, then σ → 0. In such a scenario, we in principal reach to the single layer case for the oscillators locked in their respective layers, thus the order parameter turns into the following: For λr l = γ, Eq. 17 yields r b lc = π 4 . Hence, the backward critical coupling strength is given by Forward critical coupling: From Eq. 11, the contribution from inter-layer coupling yields bounds (maximal or minimal) of ±α/2. Moreover, in the incoherent state at λ → λ f lc (r l = 0), the contribution from intra-layer mean-field is negligible. Hence, the evolution of the nodes at λ f lc is driven entirely by the effective critical frequency ω So, the critical frequencies at λ f lc are bounded in the effective interval ω lc ∈ (−γ − α/2, γ + α/2). Hence, the dynamics of the nodes at λ → λ f lc can be approximated aṡ θ Next, following the same methodology adopted for the backward transition, the forward critical threshold is given by In Fig. 6, numerical results for the order parameter and the forward and backward thresholds for either GC layers (since both the GC layers synchronizes simultaneously) and their analytical predictions given by Eqs. 16, 18 and 21, are shown for different values of α when = 0.1. The analytical predictions match fairly well with their numerical estimations. From theoretical and numerical outcomes we deduce that the threshold for first-order de-synchronization does not depend on either α or ε, however the threshold for the onset of first-order transition to synchronization does depend on α. Further, we show the dependence of λ f lc and λ b lc on factor α both numerically as well as analytically as shown in Fig. 7. The α − λ phase plot unveils that λ b lc remains independent of α and fixed to 4γ π . However, λ f lc does show dependence on α closely following Eq. 21.

Robustness against network size N and structure
We explored the phenomena of ES as a consequence of adaptive multiplexing for larger size as well as for different network architecture. Numerical results for larger sizes N = 10, 000 of multiplex networks composed of two different pairs of homogeneous topology, namely ER-ER and GC-GC network are shown in Fig. 8. Numerical results for the order parameter corresponding to different values of α for large size N are consistent with those obtained for rather small network size N = 1, 000. Hence, the occurrence of ES as an outcome of inter-layer adaptation is robust against networks size N and homogeneous network topology of the multiplex network. Nevertheless, the employed inter-layer Hebbian adaptive mechanism is incapable in triggering ES transition in heterogeneous topology for a multiplexed layer as can be observed for BA layer in the phase plots for ER-BA multiplex configuration in Fig. 5.

Conclusion
We studied the adaptive evolution of connection weights of inter-layer links in multiplex networks. The connection weight between a pair of interconnected nodes is strengthened if they are in phase and weakened if they are out of phase. Such an Hebbian learning rule plays an important role in shaping the dynamics of the phases, as well as inter-layer links control, in turn, intra-layer synchronization. In the asynchronous state of each layer, the existing Hebbian learning adaptation divides the inter-layer weight population into two almost equal groups corresponding to steady states α and −α, yielding a phase-difference of 0 and π, respectively, between the interconnected nodes. It is the almost half antiweight population −α the one which generates frustration at its end interconnected nodes and induces a first-order like transition. In the synchronous state of each layer, the Hebbian learning adaptation pulls the entire inter-layer link population into a single group corresponding to a −α steady state with the interconnected nodes maintaining a phase-difference of π. It is also unveiled that Hebbian learning rule provides a great amount of control over the abruptness (of the first-order transition) and the size of the associated hysteresis by means of the amplification factor α, however the effect of the learning rate ε is weak. The backward critical λ b c is independent on both α and ε, while λ f c is shown to have explicit dependence on α. The numerical estimation of both λ f c and λ b c have shown to fall in good agreement with their respective analytical predictions. The proposed scheme of inter-layer Hebbian adaptation is capable to bringing about first-order transition in homogeneous multiplexed layers. Phase diagrams for hysteresis width |r f l − r b l | in α − λ space for different rates ε are provided for multiplex networks with different combinations of topology.
The microscopic dynamics behind the origin of ES is simplified further when taking into account ω 2 for multiplex networks made of two ER layers. The first column in Fig. A1 shows the absence of neutral cluster around D [ii] → 0 which is present for the case of ω 2 , hence the identical frequencies of the interconnected nodes leads to only two pure anti-weight states, namely D ii α and D ii −α. And for that matter, only two equal sized clusters are obtained following |θ 1 | = 0 and π in the incoherent state (see second column). Anyhow, one single peak is obtained in the coherent state. In the third column, the behavior of the final states of D [ii] against ∆θ [i] affirms the fact that for any α in the coherent state, a single phase-cluster with |θ 1 | = π in each layer. The fourth column infers that the abrupt jump in R 2 at the onset of synchronization shows the existence of anti-phase layers. Also, the special case of ω 2 shows a handsome jump size in R 2 for any value of α as that for the case of ω  1 | belonging to the incoherent and coherent state, D [ii] against ∆θ [i] , and R 1 and R 2 order parameter as a function of λ for the multiplex network comprising two ER layers. Results are produced for N =1000 nodes, ε = 0.1 and γ = 0.5 and different values of α.