During the early eighteenth century, Leonhard Euler published a paper on The Seven Bridges of Königsberg providing a mathematical background on vertices and edges [1]. Later on, it became the cornerstone of the field of network science. Network science presents a unique platform to study various complex real-world systems by analyzing the interactions between its constituent entities and collectively investigating its behaviors [2–4]. A recent addition to the network science is the multiplex framework which incorporates multiple types of interactions among nodes by representing them in different layers [5, 6]. For example, the neurons in the brain form different groups consisting of the same neurons but interacting in different ways (chemical interaction or electrical synapses) to perform different tasks [7, 8]. Multiplex framework divides these neuronal groups into different layers based on their functionalities [9]. Similarly, transportation networks, communication networks, social networks and a lot of other real-world networks can be represented in a multiplex framework to understand their structural and dynamical features in a better fashion [10].

Recently, various synchronization scenarios have been investigated in multilayer structures, including remote and relay synchronization [11–14]. Moreover, it has been shown that multiplexing can be used to control spatio-temporal patterns in networks [15–18]. The advantage of control schemes based on multiplexing is that they allow to achieve the desired state in a certain layer without manipulating its parameters, and they can work for weak inter-layer coupling. For example, it has been found that weak multiplexing can induce coherence resonance [19] as well as chimera states [18] in neural networks.

Chimera state is a peculiar partial synchronization pattern that refers to a hybrid dynamics where coherence and incoherence emerge simultaneously in a network of identical oscillators [20–22]. Since its inception [23, 24], chimera state has attracted massive interest from the nonlinear community for both its significance in understanding complex spatiotemporal patterns and its probable applicabilities in various fields, especially in neuroscience [25]. Here we study the role of the interplay of intra- and inter-layer time delays for the emergence of chimera states in a multi(tri)plex network. Time delays represent an essential factor in real-world networks due to the finite speed of information propagating through channels connecting the nodes. They play a crucial role in determining the dynamical behavior of a complex system [26–32]. Time delays have been shown to heavily influence the parameter range for which chimera states appear for both single and multiplex networks [14, 15, 33–35]. Moreover, recently a scheme has been proposed for engineering chimera states using suitably placed heterogeneous delays [36]. In the present work we demonstrate that just the variation of delay values in the intra- and inter-layer edges can lead to various dynamical states and allows for control of spatio-temporal patterns.

We study a multiplex network with three layers (triplex) as shown in Figure 1. Every single layer represents a ring of N identical FitzHugh-Nagumo (FHN) oscillators with non-local (intra-layer) topology. The outer layers i = 1 and i = 3 are coupled with the middle layer i = 2, so that it acts as a relay layer. There is no direct connection between layers 1 and 3 (so-called ordinal coupling). Our system is described by the following equations:

where x=(u,v)T∈ℝ2, k ∈ {1, …, N}, i ∈ {1, …, 3} with all indices modulo N, describe the set of activator (u) and inhibitor (v) variables. The intra-layer delay time is τ_i and the inter-layer delay time is τ_ij. The coupling radius in layer i is given by R_i. The local dynamics of each oscillator is given by

where ε = 0.05 is the parameter characterizing the time scale separation. The FHN oscillator exhibits either oscillatory (|a| < 1) or excitable (|a| > 1) behavior depending on the threshold parameter a. In this work we focus on the oscillatory regime (a = 0.5). The parameter σ_i stands for the coupling strength inside the layer (intra-layer coupling), and σ_ij is the inter-layer coupling. For an ordinal inter-layer coupling with constant row sum we set σ₁₂ = σ₂₃, which yields the inter-layer coupling matrix

The connections between the nodes are given by the diffusive coupling with the following coupling matrix

and coupling phase ϕ=π2-0.1 [37]. This coupling configuration (i.e., predominantly activator-inhibitor cross-coupling) is similar to a phase-lag of approximately π/2 in the Kuramoto model that ensures the occurrence of chimera states [37].

Chimera states are spatio-temporal patterns where incoherent and coherent domains coexist in space. For certain values of coupling strength σ_i and coupling radius R_i one can detect them in the ith layer [37]. Recently, the phenomenon of relay synchronization of chimera states has been studied in a three-layer network of FHN oscillators [38]. In more detail, for varying the coupling delay and strength in the inter-layer connections relay synchronization of chimera states in the outer network layers has been reported. For appropriate parameters the so-called “double” chimeras are possible where the coherent parts of the chimera states are synchronized, whereas the incoherent parts remain desynchronized. Additionally, the transitions between different synchronization scenarios have been studied. Moreover, time delay in the inter-layer coupling has been shown to be a powerful tool for controlling various partial synchronization patterns in the three-layer network [14]. Therefore, by multiplexing and introducing inter-layer time delays it is possible to destroy or induce chimera states. Here we study the interplay of inter- and intra-layer time delay.

To provide an overview of the patterns observed in the network, we calculate the map of regimes in the parameter plane of intra-layer delay time τ_i and inter-layer delay time τ_ij (Figure 2). The dominating region is the one corresponding to coherent states (blue region in Figure 2). On the one hand, we detect the in-phase synchronization regime (see Figure 3C), on the other hand, we also observe a region of coherent traveling waves (see Figure 4C). By varying the delay times we can not only switch between these states, but also adjust the speed of traveling waves. In addition, we can observe salt and pepper states (green region in Figure 2), where all nodes oscillate with the same phase velocity but they are distributed between states with phase lag π incoherently [40] (see Figures 3B, 4B). The reason for this are strong variations on very short length scales, so that the dynamical patterns have arbitrarily short wavelengths. Besides these two patterns characterized by the same mean phase velocity for all the nodes in the network, we also observe chimera states (red region in Figure 2), where the oscillators within each layer show a characteristic arc-shaped mean phase velocity profile. The mean phase velocities of the oscillators are given by ω_k = 2πS_k/ΔT, k = 1, …, N, where S_k denotes the number of complete rotations performed by the kth oscillator during the time ΔT. The value ΔT = 10, 000 is fixed throughout the paper. Because of the delay, the system often exhibits traveling structures. Therefore, the classical arc-shaped profile is transformed into a wider, conical one [41] (see Figures 3A, 4A (right column)). We distinguish between the different regions on the one hand, by analyzing the mean phase velocity and a snapshot of variables u_k, on the other hand, by means of the Laplacian distance measure [39]. The latter identifies strong local curvature in an otherwise smooth spatial profile of xk by calculating the discrete Laplacian ∥(xk+1-xk)-(xk-xk-1)∥ for each k. In case of coherent dynamics we obtain low values (≈ 0) for all k, in case of salt and pepper states we get high values (> 2). Chimera states show low and high values for the coherent and incoherent domains, respectively. By taking the mean over all k, we can distinguish between coherent, chimera and “salt & pepper” states. For the numerical integration an Euler integration method is used with a step size Δt = 0.0025 and a transient time t_trans = 2, 000. All simulations are evaluated then after ΔT = 10, 000 as mentioned above.

In many systems with time delays resonance effects can be expected if the delay time is an integer or half-integer multiple of the period in the uncoupled system [42–44]. Regarding the inter-layer delay time τ_ij we can observe this effect for half-integer multiples of the period T = 2.3 of a uncoupled FHN oscillator (see Figure 2 for τ_i ≈ 2.5). Concerning the intra-layer delay time τ_i we find a resonance effect in the case of integer multiple of the delay. For greater values of the delay times τ_i and τ_ij the dynamical regions are becoming curved (see green islands in Figure 2 at τ_i ≈ 2.5). This can be explained by the fact that branches of periodic solutions, which are reappearing for integer multiples of the intrinsic period, are becoming stretched with increasing delay time [41]. In comparison to the almost vertical shape at τ_i ≈ 0.5, the green islands at τ_i ≈ 2.5 are rotated clockwise by approximate π/8. The consequence is an overlapping of the delay islands for small intra-layer delay time τ_i, whereas for larger delays the islands become separated.

Networks with multiple layers demonstrate remote synchronization of distant layers via a relay layer. Regarding the inter-layer synchronization, two synchronization mechanisms are conceivable in a triplex network:

A useful measure for synchronization between two layers i, j is given by the global inter-layer synchronization error E^ij [14, 45]:

where ∥·∥ denotes the Euclidean norm. One can distinguish between the two synchronization mechanisms by measuring the global inter-layer synchronization error between the first and second layer E¹² and between the first and third layer E¹³: In the case of full inter-layer synchronization E¹² = 0 as well as E¹³ = 0, while in the case of relay inter-layer synchronization E¹² ≠ 0 and E¹³ = 0.

To provide more insight into the synchronizability of patterns between two layers i, j, the local inter-layer synchronization error in dependence of every single node k can be used [14]:

The local inter-layer synchronization error is convenient for detecting the synchronized nodes between two layers. In Figures 3, 4 (right column) Ekij is plotted (light yellow) together with the mean phase velocity (blue): Depending on the delay times τ_i and τ_ij we can find full inter-layer synchronization (see Figure 3C) as well as relay inter-layer synchronization (see Figures 3A,B, 4A,C). These synchronization scenarios can be found for both coherent and incoherent dynamics. An additional effect is the partial relay synchronization scenario in Figure 4B: In all three layers “salt & pepper” dynamics can be observed. The nodes in the outer layers are almost all synchronized, but a small part of them destroys the relay synchronization. On the other hand a few oscillators in the relay layer are synchronized with the outer layers.

In conclusion, we have studied chimera states in a three-layer network of FitzHugh-Nagumo oscillators, where each layer has a nonlocal coupling topology. Focusing on the role of time delays in the coupling terms and their influence on chimera states, we have performed a numerical study of complex spatio-temporal patterns in the network. In the parameter plane of the intra-layer τ_i and the inter-layer τ_ij time delay, we have determined the regions where chimera patterns occur, alternating with regimes of coherent states. A proper choice of time delay allows to achieve the desired state of the network: chimera state or coherent pattern, full or relay inter-layer synchronization.

Combining the delayed interactions with the multiplex framework considered in this work can provide additional insight into the formation of the complex spatio-temporal patterns in real-world systems. Specifically, in brain networks where EEG patterns are recently reported to display chimera-like behavior at the onset of a seizure [46–48]. Inducing the chimera states by tuning the inter- and intra-layer delay values provides us with a powerful tool to control chimera states.

JS did the numerical simulations and the theoretical analysis. AZ supervised the study. All authors designed the study and contributed to the preparation of the manuscript. All the authors have read and approved the final manuscript.