Interlayer Hebbian plasticity induces first-order transition in multiplex networks

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)

Zoom In Zoom Out Reset image size

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

$$\begin{equation}\begin{aligned}\hfill {\dot {\theta }}_{1}^{\left[i\right]}& ={\omega }_{1}^{\left[i\right]}+{\lambda }_{1}\sum _{j=1}^{N}{A}_{1}^{\left[ij\right]}\enspace \mathrm{sin}\left({\theta }_{1}^{\left[j\right]}-{\theta }_{1}^{\left[i\right]}\right)+{D}^{\left[ii\right]}\enspace \mathrm{sin}\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right),\hfill \\ \hfill {\dot {\theta }}_{2}^{\left[i\right]}& ={\omega }_{2}^{\left[i\right]}+{\lambda }_{2}\sum _{j=1}^{N}{A}_{2}^{\left[ij\right]}\enspace \mathrm{sin}\left({\theta }_{2}^{\left[j\right]}-{\theta }_{2}^{\left[i\right]}\right)+{D}^{\left[ii\right]}\enspace \mathrm{sin}\left({\theta }_{1}^{\left[i\right]}-{\theta }_{2}^{\left[i\right]}\right),\hfill \end{aligned}\end{equation} \tag{ 1 }$$

where i = 1, ..., N, subscripts 1 and 2 stand for the two distinct layers, and λ represents intra-layer coupling strength, here λ₁ = λ₂ = λ. The instantaneous phase of the ith node is denoted by ${\theta }_{l}^{\left[i\right]}\left(t\right)$ and its natural frequency ${\omega }_{l}^{\left[i\right]}$ 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}_{l}^{\left[ij\right]}=1\enspace \left(0\right)$ if ith and jth 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

$$\begin{equation}{\dot {D}}^{\left[ii\right]}=\varepsilon \left[\alpha \enspace \mathrm{cos}\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right)-{D}^{\left[ii\right]}\right],\end{equation} \tag{ 2 }$$

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:

$$\begin{equation}\mathcal{M}\left(t\right)=\left(\begin{matrix}\hfill {A}_{1}\hfill & \hfill D\left(t\right)I\hfill \\ \hfill D\left(t\right)I\hfill & \hfill {A}_{2}\hfill \end{matrix}\right),\end{equation} \tag{ 3 }$$

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

$$\begin{equation}{r}_{l}\left(t\right){e}^{{\imath}{\psi }_{l}}=\frac{{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}{e}^{{\imath}{\theta }_{l}^{\left[i\right]}}}{{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}}.\end{equation} \tag{ 4 }$$

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

$$\begin{equation}{R}_{1}\left(t\right){e}^{{\imath}{\Psi}}=\frac{{\sum }_{l=1}^{2}{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}{e}^{{\imath}{\theta }_{l}^{\left[i\right]}}}{{\sum }_{l=1}^{2}{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}}.\end{equation} \tag{ 5 }$$

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

$$\begin{equation}{R}_{2}=\vert R-{R}_{1}\vert ,\end{equation} \tag{ 6 }$$

where

$$\begin{equation*}R\left(t\right){e}^{{\imath}{\Phi}}=\frac{{\sum }_{l=1}^{2}{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}{e}^{2{\imath}{\theta }_{l}^{\left[i\right]}}}{{\sum }_{l=1}^{2}{\sum }_{i=1}^{N}{k}_{l}^{\left[i\right]}}.\end{equation*}$$

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.

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₁⟩ = ⟨k₂⟩ = 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 ${\lambda }_{\mathrm{c}}^{\mathrm{b}}$ is independent of the rate ɛ. However, the forward critical coupling strength ${\lambda }_{\mathrm{c}}^{\mathrm{f}}$ 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.

Zoom In Zoom Out Reset image size

Next, 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 ${\lambda }_{\mathrm{c}}^{\mathrm{f}}$ significantly increases while the backward critical coupling ${\lambda }_{\mathrm{c}}^{\mathrm{b}}$ remains independent on α as well. Thus, each increase in α leads to a broader hysteresis width.

Zoom In Zoom Out Reset image size

To 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 ${\Delta}{\theta }^{\left[i\right]}=\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert$, global in-phase (R₁) and anti-phase (R₂) order parameter in figure 4. Note that the similar investigations are also presented for the case of ${\omega }_{1}^{\left[i\right]}={\omega }_{2}^{\left[i\right]}$ in appendix A. The first column of figure 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 D^[ii] = D^[ii](t) are bounded in the interval [−α, α]. In the incoherent state ($\lambda {< }{\lambda }_{\mathrm{c}}^{\mathrm{f},\mathrm{b}}$), 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 ($\lambda { >}{\lambda }_{\mathrm{c}}^{\mathrm{f},\mathrm{b}}$). 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 ${\lambda }_{\mathrm{c}}^{\mathrm{f}}$ with each increase in α. For transition to desynchronization, the behavior of D^[ii] when in synchronous and asynchronous state, is explained in appendix B.

Zoom In Zoom Out Reset image size

In the second column of figure 4, we study the distribution P(Δθ^[i]) of ${\Delta}{\theta }^{\left[i\right]}=\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert$, the difference between phases of the interconnected nodes for different values of α while keeping ɛ fixed. It unveils that in the incoherent state ($\lambda {< }{\lambda }_{\mathrm{c}}^{\mathrm{f}}$) belonging to an intermediate or higher value of α, two phase-clusters in each layer exist: one corresponding to ${\theta }_{2}^{\left[i\right]}={\theta }_{1}^{\left[i\right]}$ and the other to $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =\pi$. 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 ($\lambda { >}{\lambda }_{\mathrm{c}}^{\mathrm{f}}$) unveils a sharp single peaked (unimodal) distribution centered at $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =\pi$ 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 B reports the details of the behavior of P(Δθ^[i]) for both the synchronous and asynchronous state.

The third column of figure 4 illustrates R₁ (equation (5)) and R₂ (equation (6)) for different values of α. In the incoherent state one has that R₁ → 0 and R₂ → 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₁ → 0, while R₂ displays an abrupt jump at the critical coupling strength. The height of the abrupt jump for R₂ increases as α increases.

Now, since ${\dot {D}}^{\left[ii\right]}=0$ (from equation (2)) in the stationary state, one obtains (for ɛ ≠ 0)

$$\begin{equation}{D}^{\left[ii\right]}=\alpha \enspace \mathrm{cos}\left[{\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right]=\alpha \enspace \mathrm{cos}\left({\Delta}{\theta }^{\left[i\right]}\right).\end{equation} \tag{ 7 }$$

Further, it is revealed from figure 4 that for any value of α, there exists one-cluster coherent state in each layer obeying $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert \simeq \pi$, 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 $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =\mathrm{arccos}\left({\pm}1\right)=0,\pi$. 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 ${\lambda }_{\mathrm{c}}^{\mathrm{f}}$ for different values of ɛ. Figure 5 reports on how the hysteresis width $\vert {r}_{l}^{\mathrm{f}}-{r}_{l}^{\mathrm{b}}\vert$ 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., ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}={\lambda }_{l\mathrm{c}}^{\mathrm{b}}$ 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.

Zoom In Zoom Out Reset image size

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

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}_{l}^{\left[i,i\right]}=1/N$ 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

$$\begin{equation}{r}_{l}{e}^{\text{i}{\psi }_{l}}=\frac{1}{N}\sum _{j=1}^{N}{e}^{{\imath}{\theta }_{l}^{\left[j\right]}}.\end{equation} \tag{ 8 }$$

Note that the intrinsic frequencies of the nodes in both GC layers are selected from a uniform distribution ${\omega }_{l}^{\left[i\right]}=-\gamma +\frac{\gamma }{N}\left(2i-1\right)$ in the interval [−γ, γ] [44, 45]. Now equation (1) can be rewritten in the mean-field form using equation (8) as

$$\begin{equation}\begin{aligned}\hfill {\dot {\theta }}_{1}^{\left[i\right]}& ={\omega }_{1}^{\left[i\right]}+\lambda {r}_{1}\enspace \mathrm{sin}\left({\psi }_{1}-{\theta }_{1}^{\left[i\right]}\right)+{D}^{\left[ii\right]}\enspace \mathrm{sin}\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right),\hfill \\ \hfill {\dot {\theta }}_{2}^{\left[i\right]}& ={\omega }_{2}^{\left[i\right]}+\lambda {r}_{2}\enspace \mathrm{sin}\left({\psi }_{2}-{\theta }_{2}^{\left[i\right]}\right)-{D}^{\left[ii\right]}\enspace \mathrm{sin}\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right).\hfill \end{aligned}\end{equation} \tag{ 9 }$$

In the stationary state ${\dot {D}}^{\left[ii\right]}=0$, for ɛ ≠ 0, thereby equation (2) leads to

$$\begin{equation}{D}^{\left[ii\right]}=\alpha \enspace \mathrm{cos}\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right).\end{equation} \tag{ 10 }$$

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

$$\begin{equation}\begin{aligned}\hfill {\dot {\theta }}_{1}^{\left[i\right]}& ={\omega }_{1}^{\left[i\right]}+\lambda {r}_{1}\enspace \mathrm{sin}\left({\psi }_{1}-{\theta }_{1}^{\left[i\right]}\right)+\frac{\alpha }{2}\enspace \mathrm{sin}\left[2\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right)\right],\hfill \\ \hfill {\dot {\theta }}_{2}^{\left[i\right]}& ={\omega }_{2}^{\left[i\right]}+\lambda {r}_{2}\enspace \mathrm{sin}\left({\psi }_{2}-{\theta }_{2}^{\left[i\right]}\right)-\frac{\alpha }{2}\enspace \mathrm{sin}\left[2\left({\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\right)\right].\hfill \end{aligned}\end{equation} \tag{ 11 }$$

Next, we gather from numerical simulations (see figure 5) that the distribution P(Δθ^[i]) of ${\Delta}{\theta }^{\left[i\right]}=\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert$ for the two GC layers in the coherent state ($\lambda { >}{\lambda }_{\mathrm{c}}^{\mathrm{b}}$) 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 $\sigma =\frac{\alpha }{2}\enspace \mathrm{sin}\left(2{\Delta}\theta \right)$, and the model equation (9) can be rewritten as

$$\begin{equation}\begin{aligned}\hfill {\dot {\theta }}_{1}^{\left[i\right]}={\omega }_{1}^{\left[i\right]}+\sigma -\lambda {r}_{1}\enspace \mathrm{sin}\left({\theta }_{1}^{\left[i\right]}-{\psi }_{1}\right),\\ \hfill {\dot {\theta }}_{2}^{\left[i\right]}={\omega }_{2}^{\left[i\right]}-\sigma -\lambda {r}_{2}\enspace \mathrm{sin}\left({\theta }_{2}^{\left[i\right]}-{\psi }_{2}\right).\end{aligned}\end{equation} \tag{ 12 }$$

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 ${\dot {\theta }}_{1}^{\left[i\right]}={{\Omega}}_{1}$ and ${\dot {\theta }}_{2}^{\left[i\right]}={{\Omega}}_{2}$, i.e., $\vert {\omega }_{1}^{\left[i\right]}\vert -{{\Omega}}_{1}+\sigma {\leqslant}\lambda {r}_{1}$ and $\vert {\omega }_{2}^{\left[i\right]}\vert -{{\Omega}}_{2}-\sigma {\leqslant}\lambda {r}_{2}$. Hence, one obtains $\mathrm{sin}\left({\theta }_{l}^{\left[i\right]}-{\psi }_{l}\right)=\frac{{\omega }_{l}^{\left[i\right]}-{{\Omega}}_{l}{\pm}\sigma }{\lambda {r}_{l}}$, 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:

$$\begin{align*}\hfill \vert {\omega }_{1}^{\left[i\right]}\vert -{{\Omega}}_{1}+\sigma {\leqslant}\lambda {r}_{1},\\ \hfill \vert {\omega }_{2}^{\left[i\right]}\vert -{{\Omega}}_{2}-\sigma {\leqslant}\lambda {r}_{2}.\end{align*}$$

The order parameter in equation (8) for the phase-locked oscillators can be expressed as

$$\begin{equation}{r}_{l}=\frac{1}{N}\sum _{\text{CondI}}\mathrm{cos}\left({\theta }_{l}^{\left[j\right]}-{\psi }_{l}\right),\end{equation} \tag{ 13 }$$

which can be further simplified as

$$\begin{equation}\begin{aligned}\hfill {r}_{l}& =\frac{1}{N}\sum _{\text{CondI}}\sqrt{1-{\left(\frac{{\omega }_{l}^{\left[i\right]}-{{\Omega}}_{l}{\pm}\sigma }{\lambda {r}_{l}}\right)}^{2}},\hfill \\ \hfill {{\Omega}}_{l}& =\frac{{\sum }_{\text{CondI}}{\omega }_{l}^{\left[i\right]}}{{\sum }_{\text{CondI}}1}.\hfill \end{aligned}\end{equation} \tag{ 14 }$$

Now, Δθ → 0 (Δθ ≃ 0.002–0.003) for any $\lambda { >}{\lambda }_{l\mathrm{c}}^{\mathrm{b}}$ corresponding to locked state, hence σ → 0 (see figure 6).

Zoom In Zoom Out Reset image size

Backward critical coupling: in the continuum limit N → ∞, the order parameter can be rewritten in its integral form

$$\begin{equation}{r}_{l}={\int }_{\text{CondI}}\mathrm{d}{\omega }_{l}\enspace g\left({\omega }_{l}\right)\sqrt{1-{\left(\frac{{\omega }_{l}-{{\Omega}}_{l}{\pm}\sigma }{\lambda {r}_{l}}\right)}^{2}}.\end{equation} \tag{ 15 }$$

For a uniform distribution $g\left({\omega }_{l}\right)=\frac{1}{2\gamma }$ for |ω_l | < γ and since Ω_l → 0, the order parameter in equation (15) reduces to

$$\begin{equation}{r}_{l}=\frac{1}{\gamma }{\int }_{0}^{\gamma }\mathrm{d}{\omega }_{l}\sqrt{1-{\left(\frac{{\omega }_{l}{\pm}\sigma }{\lambda {r}_{l}}\right)}^{2}}.\end{equation} \tag{ 16 }$$

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:

$$\begin{equation}{r}_{l}=\frac{1}{2}\sqrt{1-{\left(\frac{\gamma }{\lambda {r}_{l}}\right)}^{2}}+\frac{\lambda {r}_{l}}{2\gamma }\enspace \mathrm{arcsin}\left(\frac{\gamma }{\lambda {r}_{l}}\right).\end{equation} \tag{ 17 }$$

For λr_l = γ, equation (17) yields ${r}_{l\mathrm{c}}^{\mathrm{b}}=\frac{\pi }{4}$. Hence, the backward critical coupling strength is given by

$$\begin{equation}{\lambda }_{l\mathrm{c}}^{\mathrm{b}}=\frac{\gamma }{{r}_{l\mathrm{c}}^{\mathrm{b}}}=\frac{4\gamma }{\pi }.\end{equation} \tag{ 18 }$$

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 $\lambda \to {\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ (r_l = 0), the contribution from intra-layer mean-field is negligible. Hence, the evolution of the nodes at ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ is driven entirely by the effective critical frequency

$$\begin{equation}{\omega }_{l\mathrm{c}}^{\left[i\right]}\simeq {\omega }_{l}^{\left[i\right]}{\pm}\alpha /2.\end{equation} \tag{ 19 }$$

So, the critical frequencies at ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ are bounded in the effective interval ${\omega }_{l\mathrm{c}}^{\left[i\right]}\in \left(-\gamma -\alpha /2,\gamma +\alpha /2\right)$. Hence, the dynamics of the nodes at $\lambda \to {\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ can be approximated as

$$\begin{equation}{\dot {\theta }}_{l}^{\left[i\right]}={\omega }_{l\mathrm{c}}^{\left[i\right]}+\lambda {r}_{l}\enspace \mathrm{sin}\left({\psi }_{l}-{\theta }_{l}^{\left[i\right]}\right).\end{equation} \tag{ 20 }$$

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

$$\begin{equation}{\lambda }_{l\mathrm{c}}^{\mathrm{f}}\simeq \frac{4}{\pi }\left(\gamma +\alpha /2\right).\end{equation} \tag{ 21 }$$

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 ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ and ${\lambda }_{l\mathrm{c}}^{\mathrm{b}}$ on factor α both numerically as well as analytically as shown in figure 7. The α − λ phase plot unveils that ${\lambda }_{l\mathrm{c}}^{\mathrm{b}}$ remains independent of α and fixed to $\frac{4\gamma }{\pi }$. However, ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ does show dependence on α closely following equation (21).

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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₁⟩ is kept fixed to 10 while ⟨k₂⟩ is tested with different values of average degree. In figure 9(a), the two ER layers with ⟨k₁⟩ = ⟨k₂⟩ = 10 synchronize and de-synchronize abruptly and simultaneously at the same critical threshold. However in figure 9(b), the layer with contrasting ⟨k₂⟩ = 20 abruptly synchronizes relatively at a lower forward threshold ${\lambda }_{2\mathrm{c}}^{\mathrm{f}}$, because of a large number of intra-layer connections, as compared to the threshold ${\lambda }_{1\mathrm{c}}^{\mathrm{f}}$ of layer with ⟨k₁⟩ = 10. Now even higher ⟨k₂⟩ = 30 sets the onset of transition at further lower threshold ${\lambda }_{2\mathrm{c}}^{\mathrm{f}}$. Moreover, abrupt de-synchronization of the two layers takes place at the same critical threshold ${\lambda }_{1\mathrm{c}}^{\mathrm{b}}={\lambda }_{2\mathrm{c}}^{\mathrm{b}}$. Not only that both ${\lambda }_{l\mathrm{c}}^{\mathrm{f}}$ and ${\lambda }_{l\mathrm{c}}^{\mathrm{b}}$ for the two layers shift towards lower values with each increasing difference in ⟨k₁⟩ and ⟨k₂⟩. Hence, contrasting average connectivities of the layers set the onset of ES transition at lower threshold.

Zoom In Zoom Out Reset image size

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.

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.

The microscopic dynamics behind the origin of ES is simplified further when taking into account ${\omega }_{1}^{\left[i\right]}={\omega }_{2}^{\left[i\right]}$ 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 ${\omega }_{1}^{\left[i\right]}\ne {\omega }_{2}^{\left[i\right]}$, 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 $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =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 $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =\pi$ 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 ${\theta }_{2}^{\left[i\right]}={\theta }_{1}^{\left[i\right]}$ in each layer while the other-half population experiencing −α leads to phase-cluster with $\vert {\theta }_{2}^{\left[i\right]}-{\theta }_{1}^{\left[i\right]}\vert =\pi$ in each layer. The fourth column infers that the abrupt jump in R₂ at the onset of synchronization shows the existence of anti-phase layers. Also, the special case of ${\omega }_{1}^{\left[i\right]}={\omega }_{2}^{\left[i\right]}$ shows a handsome jump size in R₂ for any value of α as that for the case of ${\omega }_{1}^{\left[i\right]}\ne {\omega }_{2}^{\left[i\right]}$.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size