## Abstract

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 interlayer 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.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

## 1. Introduction

The study of synchronization, a collective motion of initially nonidentical 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 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 (ES)) 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 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–10] and noise [10, 11], or by the presence of a fraction of adaptively coupled oscillators [12–17], mean-field diffusion [18], traffic processes [19], the presence of nearest-neighbor competitive interaction or symmetry-breaking 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–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–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 segregates the population of interlayer weights into almost equal number of excitatory and inhibitory weights for a 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. (figure 1)

## 2. 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 and its natural frequency 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 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 rule

where *α* ∈ [0, 1] is a factor which amplifies the amount of learning if the two nodes are synchronized and *ɛ* ∈ [0, 1] is the learning rate. The *D*^{[ii]} in equation (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 order parameter *r*_{l
} for a layer *l* in terms of the average phase *ψ*_{l
} as

*r*_{l
} = 1 represents a completely synchronous state, while *r*_{l
} = 0 implies total incoherence. In a similar way, one can define global in-phase (one-cluster) order parameter for the entire multiplex network as

Furthermore, to capture the degree of global anti-phase (two-cluster) synchronization, a new global order parameter *R*_{2} (dipole moment of the distribution of anti-phases) [35, 41] is defined as follows

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 equation (6) for one-cluster synchronization.

## 3. Results

We consider a multiplex network made of two different realizations of 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 different realizations of initial phases and natural frequencies drawn from a uniform random distribution such that *θ*^{[i]} ∈ [0, 2*π*) and *ω*^{[i]} ∈ [−*γ*, *γ*) where *γ* = 0.5 unless otherwise stated, respectively. The initial *N* values of *D*^{[ii]}(*t*) are selected as *D*^{[ii]}(*t* = 0) = *d*^{[i]}/*N*, where *d*^{[i]} 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. Figure 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 is independent of the rate *ɛ*. However, the forward critical coupling strength 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 hysteresis.

Download figure:

Standard image High-resolution imageNext, the impact of the amplification factor *α* with a fixed choice of *ɛ* is reported in figure 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 significantly increases while the backward critical coupling remains independent on *α* as well. Thus, each increase in *α* leads to a broader hysteresis width.

Download figure:

Standard image High-resolution imageTo gather information at a microscopic level on the origin of the first-order transition, we investigate the behavior of *D*^{[ii]}, the distribution *P*(Δ*θ*^{[i]}) of , global in-phase (*R*_{1}) and anti-phase (*R*_{2}) order parameter in figure 4. Note that the similar investigations are also presented for the case of in appendix *D*^{[ii]} in the forward continuation of *λ* for different values of *α* at a fixed rate *ɛ*. It is apparent that the stationary values of *D*^{[ii]} = *D*^{[ii]}(*t*) are bounded in the interval [−*α*, *α*]. In the incoherent state (), 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 (). 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 with each increase in *α*. For transition to desynchronization, the behavior of *D*^{[ii]} when in synchronous and asynchronous state, is explained in appendix

Download figure:

Standard image High-resolution imageIn the second column of figure 4, we study the distribution *P*(Δ*θ*^{[i]}) of , the difference between phases of the interconnected nodes for different values of *α* while keeping *ɛ* fixed. It unveils that in the incoherent state () belonging to an intermediate or higher value of *α*, two phase-clusters in each layer exist: one corresponding to 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 () unveils a sharp single peaked (unimodal) distribution centered at 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. As for the desynchronization transition, appendix *P*(Δ*θ*^{[i]}) for both the synchronous and asynchronous state.

The third column of figure 4 illustrates *R*_{1} (equation (5)) and *R*_{2} (equation (6)) for different values of *α*. In the incoherent state one has that *R*_{1} → 0 and *R*_{2} → 0, as there exists two anti-phase 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 has *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 (from equation (2)) in the stationary state, one obtains (for *ɛ* ≠ 0)

Further, it is revealed from figure 4 that for any value of *α*, there exists one-cluster coherent state in each layer obeying , hence *D*^{[ii]} ≃ −*α* from equation (7). Also, there exists two populations of inter-layer anti-weights *D*^{[ii]} ≃ *α* and *D*^{[ii]} ≃ −*α* in the incoherent state, hence equation (7) yields . 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 figure 2 highlight the 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 for different values of *ɛ*. Figure 5 reports on how the hysteresis width 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. There are systems in which the first-order transition arises without hysteresis, i.e., as well [21, 44]. Nevertheless, in our case the first-order transition arises exhibiting a good hysteresis width complementing the abrupt jump, which gradually decreases in size as *α* reduces from 1 to a lower value. A further decrement in *α* gradually switches the route of transition from a discontinuous to a continuous one with a decrease in the abrupt jump size as well as the hysteresis width. For a better subtlety, the phase plots in figure 5 are simulated for an increment and a decrement in *λ* of Δ*λ* = 0.002, and an increment in *α* of Δ*α* = 0.02.

Download figure:

Standard image High-resolution imageThe 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 because of the heterogeneity in degrees. For the same reason, a BA-BA configuration of multiplexes features a continuous route to synchronization for any value of *α* and *ɛ*.

## 4. 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 in model equation (1) and the order parameter defined in equation (4) for a complex network can then be expressed for a GC layer as following

Note that the intrinsic frequencies of the nodes in both GC layers are selected from a uniform distribution in the interval [−*γ*, *γ*] [44, 45]. Now equation (1) can be rewritten in the mean-field form using equation (8) as

In the stationary state , for *ɛ* ≠ 0, thereby equation (2) leads to

Hence, evolution of the nodes in the stationary state is ruled by the following self-consistent equations

Next, we gather from numerical simulations (see figure 5) that the distribution *P*(Δ*θ*^{[i]}) of for the two GC layers in the coherent state () 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 , and the model equation (9) can be rewritten as

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 and , i.e., and . Hence, one obtains , 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 equation (8) for the phase-locked oscillators can be expressed as

which can be further simplified as

Now, Δ*θ* → 0 (Δ*θ* ≃ 0.002–0.003) for any corresponding to locked state, hence *σ* → 0 (see figure 6).

Download figure:

Standard image High-resolution image
*Backward critical coupling:* in the continuum limit *N* → **∞**, the order parameter can be rewritten in its integral form

For a uniform distribution for |*ω*_{l
}| < *γ* and since Ω_{l
} → 0, the order parameter in equation (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
} = *γ*, equation (17) yields . Hence, the backward critical coupling strength is given by

*Forward critical coupling:* from equation (11), the contribution from inter-layer coupling yields bounds (maximal or minimal) of ±*α*/2. Moreover, in the incoherent state at (*r*_{l
} = 0), the contribution from intra-layer mean-field is negligible. Hence, the evolution of the nodes at is driven entirely by the effective critical frequency

So, the critical frequencies at are bounded in the effective interval . Hence, the dynamics of the nodes at can be approximated as

Next, following the same methodology adopted for the backward transition, the forward critical threshold is given by

In figure 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 equations (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 and on factor *α* both numerically as well as analytically as shown in figure 7. The *α* − *λ* phase plot unveils that remains independent of *α* and fixed to . However, does show dependence on *α* closely following equation (21).

Download figure:

Standard image High-resolution image## 5. 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 figure 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 figure 5.

Download figure:

Standard image High-resolution image## 6. Effect of different average connectivity ⟨*k*⟩

So far we have examined the behavior of phase transition for the same average connectivity for the two multiplexed layers. Here we examine the nature of phase transition for multiplex networks with contrasting average connectivity for the two ER layers. Figure 9 exhibits the behavior of order parameter when ⟨*k*_{1}⟩ is kept fixed to 10 while ⟨*k*_{2}⟩ is tested with different values of average degree. In figure 9(a), the two ER layers with ⟨*k*_{1}⟩ = ⟨*k*_{2}⟩ = 10 synchronize and de-synchronize abruptly and simultaneously at the same critical threshold. However in figure 9(b), the layer with contrasting ⟨*k*_{2}⟩ = 20 abruptly synchronizes relatively at a lower forward threshold , because of a large number of intra-layer connections, as compared to the threshold of layer with ⟨*k*_{1}⟩ = 10. Now even higher ⟨*k*_{2}⟩ = 30 sets the onset of transition at further lower threshold . Moreover, abrupt de-synchronization of the two layers takes place at the same critical threshold . Not only that both and for the two layers shift towards lower values with each increasing difference in ⟨*k*_{1}⟩ and ⟨*k*_{2}⟩. Hence, contrasting average connectivities of the layers set the onset of ES transition at lower threshold.

Download figure:

Standard image High-resolution image## 7. 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. Such an Hebbian learning rule plays an important role in shaping the dynamics of the phases, as well as interlayer links control, in turn, intra-layer synchronization. In the asynchronous state of both the layers, the existing Hebbian learning adaptation segregates the interlayer link population into two almost equal groups of excitatory and inhibitory weights. The half population of inhibitory (excitatory) weights makes their respective end mirror nodes anti-phase (in-phase) in the two layers. It is for the half inhibitory weight population, which causes the frustration among the nodes in the two layers and leads to a first-order transition route to synchronization. In the synchronous state of the two layers, the Hebbian learning adaptation drives the entire interlayer link population to form a global cluster of inhibitory weights, which yields two anti-phase layers. Thus the interlayer Hebbian plasticity not only triggers the first-order transition but also provides a great amount of control over the abruptness and area of the associated hysteresis loop by means of the Hebbian learning parameters. The same has been conveyed through the phase diagrams presenting hysteresis width for multiplex networks with different combinations of topology. It has also been demonstrated theoretically as well as numerically that the backward critical coupling strength is independent of both the learning parameters, while the forward critical coupling strength is shown to have explicit dependence on the learning parameters. However, the proposed Hebbian plasticity based scheme is capable of bringing about the first-order transition in only homogeneous multiplexed layers.

## Acknowledgments

SJ thanks CSIR grant 25(0293)/18/EMR-II for financial support. ADK acknowledges Govt of India CSIR grant 25(0293)/18/EMR-II for RA fellowship.

## Appendix A.: case

The microscopic dynamics behind the origin of ES is simplified further when taking into account for multiplex networks made of two ER layers. The first column in figure A1 shows the absence of neutral cluster around *D*^{[ii]} → 0 which is present for the case of , 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 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 originates from the cloud of *D*^{[ii]} → −*α*. In the incoherent state, one-half population of the inter-layer link experiencing *D*^{[ii]} → *α* gives birth to a phase-cluster with in each layer while the other-half population experiencing −*α* leads to phase-cluster with 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 shows a handsome jump size in *R*_{2} for any value of *α* as that for the case of .

Download figure:

Standard image High-resolution image## Appendix B.: Weight and phase distribution in backward continuation

Figure B1 illustrates the behavior of stationary values of *D*^{[ii]} as a function of *λ*, and the distribution *P*(Δ*θ*^{[i]}) for different values of *α* in the backward continuation. Beginning with a synchronous state the distribution *P*(Δ*θ*^{[i]}) sports a unimodal distribution centered at *π*, and even after abrupt desynchronization it tends to follow unimodal distribution with its mean still at *π* but sporting relatively a shorter peak and broader deviation. In both the asynchronous and synchronous state, a large population of interlayer links has *D*^{[ii]} → −*α*. However the fraction of interlayer links favoring *D*^{[ii]} → *α* is less and varies depending on the value of *λ*. It is approximately less than 10% and 1% of number of interlayer links in the asynchronous and synchronous state, respectively.

Download figure:

Standard image High-resolution image