Dragon Intermittency at the Transition to Synchronization in Coupled Rulkov Neurons

1. Introduction

Chaotic synchronization is a fundamental phenomenon in dynamical systems that governs the collective behavior of interacting nonlinear oscillators [1,2]. It covers a wide range of natural and engineered processes, from biological rhythms to advanced communication systems. Among the various facets of this phenomenon, transitions from asynchronous to synchronous dynamics in coupled systems are of particular importance [3,4,5,6,7]. These transitions often reveal rich and complex dynamical behaviors, including intermittent states when the system dynamics switches between order and chaos. The intermittency route to chaos is an important research topic in nonlinear dynamics, describing how systems transition to chaotic behavior as control parameters cross critical thresholds. Characterized by the alternation between nearly regular “laminar” phases and irregular “turbulent” bursts, this transition is a hallmark of many nonlinear systems. The turbulent bursts typically exhibit durations weakly dependent on the control parameter, whereas the mean duration of laminar phases decreases as the system progresses further into chaos, eventually vanishing altogether. Traditionally, such transitions are governed by singular bifurcation mechanisms that delineate the onset of chaos.

Classical forms of intermittency, as categorized by Pomeau and Manneville [8,9], provide a theoretical foundation for these transitions. Type-I intermittency emerges from saddle-node bifurcations, type-II intermittency from Hopf bifurcations, and type-III intermittency from inverse period-doubling bifurcations. Other types, such as crisis-induced intermittency, occur due to the collision of a chaotic attractor with a periodic orbit or its basin boundary [10], while multistate intermittency is induced by noise in a system with coexisting attractors [11,12]. Each type exhibits distinct scaling laws near the critical point, offering insight into the universal features of chaotic systems. For example, on-off intermittency—a common phenomenon near synchronization boundaries—displays a mean laminar phase duration that scales with the control parameter as $-1$ and a probability density of laminar durations that scales as $-3/2$ [13,14].

In coupled chaotic systems, intermittency plays a pivotal role in the transitions between asynchronous and synchronous states. Over the last two decades, extensive research has been dedicated to uncovering intermittent behaviors near synchronization regimes [15,16,17]. These studies highlight the sensitivity of intermittency to various factors, including the type of synchronization, parameter mismatches, noise, and coupling mechanisms. Near synchronization boundaries, systems frequently exhibit phenomena such as type-I intermittency [9,18,19], eyelet intermittency [20,21], or ring intermittency [22,23]. Eyelet intermittency typically occurs near the boundary separating synchronization and type-I intermittency, while ring intermittency arises in systems with significant parameter mismatches, leading to detuning of fundamental frequencies.

Initially, it was believed that coupled chaotic oscillators primarily transitioned to synchronization through mechanisms like on-off intermittency [24] or jump intermittency [25]. Jump intermittency, characterized by its irregular dynamics and intricate statistical properties, is particularly prominent near generalized synchronization. However, ongoing investigations continue to uncover new intermittency mechanisms in both continuous and discrete-time systems, reflecting the diverse pathways by which coupled chaotic systems approach synchronization. This rich tapestry of intermittency behaviors underscores the complexity of transitions in nonlinear systems and highlights the need for further exploration to unravel the underlying principles.

In this paper, we focus on the dynamics of two mutually coupled, non-identical Rulkov neurons, investigating the emergence of a novel type of intermittency during transitions to chaos. The Rulkov map, a discrete-time nonlinear model, is widely recognized for its ability to reproduce spiking and bursting dynamics observed in biological neurons with remarkable computational simplicity [26]. This makes the Rulkov model particularly suitable for exploring complex neural behaviors, including synchronization, chaos, and multistability, under varying coupling strengths and parameter mismatches. Recent studies have demonstrated that coupled Rulkov neurons exhibit intricate dynamics, such as multistability, periodicity, and chaotic transitions, depending on the control parameters [12,27].

Here, we extend this understanding by uncovering a unique intermittency phenomenon that arises when the coupled neurons attempt to synchronize but are disrupted by desynchronizing dynamics due to parameter mismatches. Unlike classical forms of intermittency, this new type—termed dragon intermittency—is distinguished by specific scaling behaviors and a probability density function with a unique "dragon-like" shape. We identify two critical power-law scalings near the onset of intermittency: (i) the mean duration of synchronization intervals scales with a control parameter with a critical exponent of $-0.5$, and (ii) the probability distribution of synchronization intervals scales with an exponent of $-1.7$ for short intervals. Interestingly, we also observe the existence of a most probable super-long laminar phase duration for each fixed coupling strength and parameter mismatch, which diminishes as the parameters are further increased.

The discovery of dragon intermittency broadens the landscape of known intermittency types in coupled dynamical systems and reveals new pathways for the interplay between synchronization and chaos in discrete-time neural models. These findings contribute to a deeper understanding of the mechanisms underlying synchronization transitions in systems with parameter mismatches, offering potential implications for both theoretical studies and practical applications, such as neural network modeling, brain dynamics, and coupled oscillator systems in physics and biology.

The remainder of this paper is organized as follows. Section 2 introduces the mathematical model of coupled Rulkov neurons and provides an overview of the numerical methods employed. Section 3 explores the observed dynamics in detail, including the coexistence of attractors and the structure of their basins of attraction. In Section 4, we present the key results, focusing on the scaling laws and statistical properties associated with dragon intermittency. Finally, Section 5 summarizes the main conclusions and suggests promising directions for future research.

2. Model of Coupled Rulkov Neurons

Consider a system of two non-identical neurons mutually coupled by electrical synapses:

$$\begin{array}{c}{x}_{t+1}=f({\gamma }_1,x_t)+\sigma (y_t-x_t),y_{t+1}=f({\gamma }_2,y_t)+\sigma (x_t-y_t),\hfill \end{array}$$

(1)

where $\sigma$ represents the coupling strength, and f is a nonlinear Rulkov function given by [26]:

$$f(\gamma ,v)={\displaystyle \frac{\alpha }{1+v^2}}+\gamma$$

with parameter $\gamma ={\gamma }_1,{\gamma }_2$, variable $v=x,y$, and constant $\alpha =4.1$. In the system (1), one of the map parameters slightly differs, such that ${\gamma }_2={\gamma }_1+\Delta$, where $\Delta$ denotes the parameter mismatch or detuning.

In the subsequent analysis, we fix ${\gamma }_1=-1.75$, a value for which the isolated map $x_{t+1}=f({\gamma }_1,x_t)$ exhibits a stable 3-cycle, as shown in the bifurcation diagram in Figure 1. For this study, we consider a small mismatch in the range $0\le \Delta \le 0.002$, ensuring that the isolated neuron is in a 3-cycle oscillatory mode.

The dynamics of the coupled system (1) are strongly influenced by the coupling parameter $\sigma$. This dependence is illustrated in Figure 2 through bifurcation diagrams, where $\sigma$ is used as the control parameter. The diagrams compare the system’s behavior in two scenarios: one with no mismatch ($\Delta =0$) and the other with slight detuning ($\Delta =0.002$). In these diagrams, the x-coordinates of various attractors, found by randomly varying initial conditions, are represented in red ($\mathcal{R}$), green ($\mathcal{G}$), and blue ($\mathcal{B}$).

In the case of an isolated one-dimensional subsystem, the phase space is one-dimensional and hosts a single attractor corresponding to a 3-cycle, as shown in Figure 1. However, in the coupled system, even when the coupling strength is zero ($\sigma =0$), the phase space becomes two-dimensional, representing two uncoupled one-dimensional subsystems. This configuration results in the coexistence of three distinct 3-cycles within the two-dimensional phase space, each corresponding to a unique synchronization mode. This multirhythmicity is clearly seen in the bifurcation diagram in Figure 2 by the $\mathcal{R}$, $\mathcal{G}$, and $\mathcal{B}$ branches.

The insets in Figure 2 provide a detailed view of these cycles, emphasizing the intricate and diverse dynamics introduced by coupling. These findings highlight the profound impact of coupling and detuning on the system’s behavior, revealing complex synchronization patterns and transitions.

At a critical value of $\sigma$, the multirhythmicity vanishes, leaving the 3-cycle $\mathcal{R}$ as the sole attractor. As $\sigma$ increases further, the system (1) transitions from the regular 3-cycle to a chaotic regime through a crisis bifurcation at ${\sigma }_{*}$. For $\Delta =0$, chaos emerges at ${\sigma }_{*}=0.0711516$ (Figure 2(a)), while for $\Delta =0.002$, it appears at ${\sigma }_{*}=0.019011$ (Figure 2(b)). Notably, this $\sigma$ threshold is significantly lower in the presence of mismatch ($\Delta =0.002$) compared to the case without mismatch ($\Delta =0$).

The results of our detailed parametric analysis of the order-to-chaos transition are presented in Figure 3, where the largest Lyapunov exponents $\Lambda$ are plotted as functions of $\sigma$ and $\Delta$. Positive values of $\Lambda$ indicate chaotic behavior. As shown in Figure 3(a), the crisis bifurcation point ${\sigma }_{*}$ decreases as the mismatch $\Delta$ is increased. In the two-parameter space, the crisis bifurcation line appears approximately linear, as depicted in Figure 3(b). These findings underscore the sensitivity of the system’s dynamics to parameter mismatches and provide valuable insights into the mechanisms driving the transition to chaos.

3. Basins of Attraction

Examine the arrangement of regular and chaotic attractors, as well as the geometry of their basins of attraction, in the phase plane. Representative examples are presented in Figure 4 and Figure 5, corresponding to $\Delta =0$ and $\Delta =0.002$, respectively.

The system (1) with zero coupling ($\sigma =0$) and no mismatch ($\Delta =0$) exhibits three stable, coexisting 3-cycles ($\mathcal{R}$, $\mathcal{G}$, and $\mathcal{B}$), as shown in Figure 4(a). The basins of attraction for these cycles form a "parquet" of rectangles: white for $\mathcal{R}$, light green for $\mathcal{G}$, and light blue for $\mathcal{B}$. When the coupling strength is slightly increased to $\sigma =0.01$ (Figure 4(b)), the attractors remain the same, but their basin structures change qualitatively, transitioning to fractal boundaries.

As $\sigma$ increases further, the stable cycles $\mathcal{G}$ and $\mathcal{B}$ vanish, and the system becomes nearly monostable with the $\mathcal{R}$ attractor dominating, as seen in the bifurcation diagram in Figure 2(a). However, narrow $\sigma$-intervals still exist where the stable $\mathcal{R}$ 3-cycle coexists with other attractors. For example, at $\sigma =0.035$ (Figure 4(c)), the $\mathcal{R}$ 3-cycle coexists with a 74-cycle, the basin of which is shown in light blue.

When $\sigma$ exceeds the crisis point ${\sigma }_{*}=0.0711516$, the stability of the $\mathcal{R}$ cycle is lost, and the system transitions to chaos. The chaotic attractor for $\sigma =0.1$ and $\Delta =0$ is shown in Figure 4(d), with a largest Lyapunov exponent of $\Lambda =0.248$.

Next, consider the effects of increasing the mismatch to $\Delta =0.002$ (Figure 5). Similar to the case with $\Delta =0$, distinct basin boundaries are observed at zero coupling ($\sigma =0$) (Figure 5(a)), and the basins exhibit fractal structures as the coupling strength increases (Figure 5(b)). However, unlike the case with $\Delta =0$, when $\Delta =0.002$, a further increase in $\sigma$ results in the annihilation of only the $\mathcal{B}$ cycle. This leaves the system bistable with two coexisting 3-cycles ($\mathcal{R}$ and $\mathcal{G}$), as shown in Figure 5(c) for $\sigma =0.01$. As $\sigma$ is increased further, the system becomes monostable with only the $\mathcal{R}$ 3-cycle for $\sigma =0.015$ (Figure 5(c)), and eventually chaotic for $\sigma >0.02$. Examples of chaotic attractors with Lyapunov exponents of $\Lambda =0.144$ (Figure 5(e)) and $\Lambda =0.267$ (Figure 5(f)) are shown for $\sigma =0.04$ and $\sigma =0.1$, respectively. As illustrated, the chaotic attractor becomes increasingly blurred as the coupling strength is raised.

Let us draw your attention to Figure 5(d), which illustrates the case of monostability. Although the attractor appears simple, it embodies a complexity not immediately apparent in this figure. While the system (1) is monostable and all initial conditions eventually converge to the same 3-cycle, $\mathcal{R}$, the convergence process itself reveals an important characteristic: the system reaches the stable solution through a chaotic transient behavior. This transient behavior is further illustrated in Figure 6, where an ensemble of solutions is shown, starting from the nodes of a uniform grid in the region $[-2,3]\times [-2,3]$. As seen in the figure, for a sufficiently long transient period, the state distribution of this ensemble is metastable and resembles the chaotic attractor shown in Figure 5(e). It takes several thousand iterations for the ensemble to begin localizing near the points of the stable 3-cycle, $\mathcal{R}$. This phenomenon, known as transient chaos [28], highlights the intricate dynamics underlying the system’s apparent simplicity.

4. Synchronization and Intermittency

We now turn our attention to how the changes in the oscillatory activity of system (1) influence synchronization dynamics. The synchronization of identical neurons ($\Delta =0$) for various values of $\sigma$ is illustrated in Figure 7 through time series. For weak coupling ($\sigma =0.01$), we observe three distinct synchronization patterns, as shown in Figure 7(a). These include complete synchronization (top panel), which occurs for initial conditions within the basin of the 3-cycle $\mathcal{R}$, as well as two forms of lag synchronization: one for solutions starting in the basins of $\mathcal{G}$ (middle panel) and another for solutions starting in the basins of $\mathcal{B}$ (bottom panel), both of which exhibit a one-step lag in synchronization.

Figure 7(b) displays two types of synchronization for $\sigma =0.035$: complete synchronization for the 3-cycle $\mathcal{R}$ (top panel) and lag synchronization for solutions starting in the basin of attraction of the 74-cycle (bottom panel). In the latter case, the synchronous dynamics are disrupted by phase slips in the x coordinate, which alternately lead it one step ahead and one step behind the y coordinate.

As $\sigma$ exceeds the critical value ${\sigma }_{*}=0.0711516$, system (1) becomes monostable with the 3-cycle $\mathcal{R}$ as the sole attractor. In this regime, the system reaches complete synchronization regardless of the initial conditions.

Synchronization in the chaotic regime is illustrated in Figure 7(c) for $\sigma =0.1$. Here, alternating time intervals of in-phase synchronization and complex mixed oscillations are observed. To identify the in-phase synchronization intervals, we introduce a synchronization indicator variable z defined as:

$$z_t=\mathrm{sgn}\left((x_{t+1}-x_t)(y_{t+1}-y_t)\right),$$

where $z=1$ indicates in-phase synchronization and $z=-1$ means asynchronous behavior. The time series of the synchronization indicator z is shown in the bottom panel of Figure 7(c). It is evident that long intervals of in-phase synchronization ($z=1$) are interrupted by alternating periods of mixed synchronous and asynchronous behaviors. One can see that even with identical neurons, increasing the coupling strength induces intermittent synchronization.

Now, let us explore the synchronization scenarios for the case of non-identical neurons with a mismatch of $\Delta =0.002$. The time series illustrating different synchronization regimes are presented in Figure 8 for various values of $\sigma$.

Figure 8(a) illustrates three variants of synchronization in the tri-rhythmic system (1) with $\sigma =0.003$, where the dynamics exhibits the coexistence of three 3-cycles ($\mathcal{R}$, $\mathcal{G}$, and $\mathcal{B}$). The basins of attraction of these cycles are shown in Figure 5(b). For stronger coupling ($\sigma =0.01$), as shown in Figure 8(b), the system transitions to a birhythmic state, where only two 3-cycles ($\mathcal{R}$ and $\mathcal{G}$) coexist, with their basins of attraction shown in Figure 5(c).

At larger coupling strengths, intermittent synchronization emerges. Figure 8(c) and Figure 8(d) display examples of this behavior, where time intervals of in-phase synchronization ($z=1$) are interrupted by periods of mixed, complex oscillations for $\sigma =0.04$ and $\sigma =0.1$, respectively. Notably, for $\sigma =0.04$, the time intervals of co-directed behavior ($\tau$) are longer than for $\sigma =0.1$.

Due to the unpredictable nature of $\tau$ in the chaotic regime, we focus on studying the probability P that the system remains in a state of co-directed oscillations ($z=1$). Figure 8(e) shows how this probability P depends on the coupling strength $\sigma$ for three different mismatch values $\Delta$. The smaller the mismatch $\Delta$, the wider the range of $\sigma$ values for which co-directed behavior ($z=1$) is observed with $P=1$. As $\sigma$ is increased, the probability P decreases. However, it is important to note that for nearly all values of $\sigma$, the probability remains greater than 0.5, indicating that in the complex chaotic regime, the time intervals of in-phase synchronization prevail over asynchronous behavior.

Let us now examine the statistical properties of the intermittency observed. First, we will analyze how the mean values $\langle \tau \rangle$ of the synchronous intervals $\tau$ vary with the parameter $\sigma$ near the crisis bifurcation point ${\sigma }_{*}$, which marks the onset of intermittency. To identify this bifurcation, Figure 9(a) and Figure 9(b) present the largest Lyapunov exponents $\Lambda$ as a function of $\sigma$ for $\Delta =0.001$ and $\Delta =0.002$, respectively. The bifurcation points ${\sigma }_{*}$ are indicated by vertical dashed lines.

Let d denotes a small deviation of the coupling parameter $\sigma$ from the crisis bifurcation point ${\sigma }_{*}$, such that $\sigma ={\sigma }_{*}+d$. Our analysis reveals that the mean duration of synchronization intervals follows a power law:

$$\langle \tau \rangle \sim d^{-0.5}.$$

In Figure 9(c) and Figure 9(d), this scaling relation is depicted by red lines for $\Delta =0.001$ and $\Delta =0.002$, respectively. The values of $\langle \tau \rangle$ obtained numerically are marked with asterisks.

Another key characteristic of intermittency is the probability distribution of time intervals $\tau$ exhibiting co-directed behavior ($z=1$). Figure 10(a) and Figure 10(b) display the time series of variables $x,y$, and synchronization indicator z for two values of deviation ($d={10}^{-4}$ and $d={10}^{-3}$) from the onset of intermittency. As shown, the time intervals of synchronous behavior become shorter and occur more frequently as d is increased.

The probability density functions $p\left(\tau \right)$ of these intervals are presented in Figure 11 for two different values of $\Delta$ and d. For small $\tau$, the asymptotic behavior follows a power law

$$p\left(\tau \right)\sim {\tau }^{-1.7},$$

which is highlighted in green.

As $\tau$ is increased, an interesting feature emerges in the probability density function $p\left(\tau \right)$. Narrow peaks with high probability become clearly visible, highlighting the most frequently observed time intervals with co-directed behavior. Notably, during these intervals, the system tends to remain near the states of the unstable 3-cycle.

5. Conclusions

In this paper, we investigated the complex dynamics of two non-identical neurons coupled through electrical synapses within a discrete dynamical framework based on the Rulkov model. Our study uncovered a rich dynamics of joint behaviors and transitions, illustrating how coupling strength and parameter mismatches govern the stability, synchronization, and complexity of the system.

We found that for small coupling strengths and parameter mismatches, the system exhibits the coexistence of three distinct 3-cycles. This tri-rhythmicity diminishes as coupling strength or mismatch parameters are increased, ultimately driving the system through a transition to chaos. This transition underscores the system’s sensitivity to both intrinsic (parameter mismatches) and extrinsic (coupling) control variables, highlighting the interplay between heterogeneity and connectivity in shaping the system dynamics.

In the chaotic regime, we identified a novel type of intermittency characterized by the statistical properties of synchronization intervals during co-directional neuron activity. Two key scaling laws emerged: (i) a power-law scaling of $\langle \tau \rangle \sim d^{-0.5}$ for the mean duration of synchronous phases near the crisis bifurcation, reflecting the system critical sensitivity near the transition to chaos and (ii) a power-law scaling of $p\left(\tau \right)\sim {\tau }^{-1.7}$, for the probability density of synchronous interval durations, revealing the statistical structure of intermittency for short intervals.

These results emphasize the presence of unique coherence properties within the intermittency dynamics. Notably, longer synchronization intervals exhibit significant coherence, a distinguishing feature of this system. The probability density function of synchronous durations also shows a distinctive "dragon-like" tail, highlighting the irregularity and unpredictability of the chaotic dynamics. This signature inspired the nomenclature for this type of intermittency and positions it as a novel phenomenon in the study of coupled systems.

The scaling laws and dynamical patterns revealed in this work enhance our understanding of fundamental mechanisms in coupled non-identical chaotic maps. Their robustness across parameter ranges suggests that these findings may generalize to other coupled chaotic systems, including neuronal networks and biological systems with similar coupling dynamics.

Future research could extend these insights by investigating the role of more complex network topologies and interactions, assessing the biological relevance of these dynamics in real neural systems, and exploring the effects of noise, time delays, or alternative coupling forms on synchronization and intermittency. By delving deeper into these universal scaling behaviors and the intricate interplay of synchronization and chaos, this study contributes to the broader understanding of coupled dynamical systems, with implications for both theoretical and applied contexts.