Stochastic Disruption of Synchronization Patterns in Coupled Nonidentical Neurons

1. Introduction

In neuronal networks, synchronization plays a crucial role in information processing, neural coding, and motor control, and its disruption has been implicated in various neurological disorders such as epilepsy, Parkinson’s disease, and schizophrenia [1,2,3,4]. Real-world neuronal networks are subject to stochastic fluctuations arising from intrinsic neural noise, synaptic variability, and external perturbations [5,6,7,8]. These fluctuations can significantly impact neuronal dynamics, leading to intermittent synchronization, chaotic transitions, and complete desynchronization [9,10,11]. In multistable systems, noise can lead to switching between coexisting attractors resulting in multistate intermittency and monostability [12,13,14,15]. Therefore, studying the stability and stochastic disruption of synchronization in coupled neurons is of great theoretical and practical significance.

In this work, we are interested in how random perturbations of the coupling parameter affect synchronization patterns in a coupled neural system. As a neural model, we consider a system of two non-identical Rulkov neurons coupled by an electrical synapse [16]. The Rulkov model is widely used to describe the bursting dynamics of neurons due to its computational efficiency and ability to capture key aspects of neuronal excitability. Various phenomena in the isolated deterministic and stochastic Rulkov model are actively studied (see, e.g. [17,18,19,20]).

Coupled Rulkov neurons exhibit rich dynamical behaviors, including multiple synchronization regimes that depend on the coupling strength and system parameters [21,22,23,24,25]. Depending on the coupling parameter, different synchronization states, including mono-, bi-, and tristable regimes, can emerge, corresponding to distinct synchronization patterns such as in-phase and lag synchronization. Through systematic parametric analysis, we identify key conditions under which noise intensity influences the transitions between different dynamical states. To analyze these complex dynamics, we employ a combination of direct numerical simulations and an analytical method based on confidence domains, which provide a robust statistical tool for assessing the variability of system trajectories under stochastic perturbations, allowing us to characterize the range of noise intensity that sustains intermittency and determine its dependence on coupling strength. The size and shape of the confidence domains are determined from the stochastic sensitivity matrices of the corresponding deterministic attractors. The stochastic sensitivity function technique was elaborated for different types of regular and chaotic attractors [26] and is now actively used in the parametric analysis of various stochastic phenomena (see, e.g. [23,27,28]).

A particularly intriguing phenomenon observed in our study is the existence of an intermittency regime, where a synchronous 3-cycle regime is interrupted by asynchronous chaotic bursts. This phenomenon underscores the critical role of basins of chaotic transients in shaping the response of the system to stochastic perturbations. The presence of intermittent chaotic behavior suggests that noise can induce transitions between metastable states, highlighting the nontrivial interplay between deterministic and stochastic influences in coupled neuronal dynamics. We investigate the stochastic disruption of lag synchronization and the emergence of a monostable in-phase synchronization regime.

The rest of this paper is organized as follows: Section 2 introduces the deterministic model of coupled Rulkov neurons and explores their dynamical behavior. Section 3 extends the discussion to the stochastic model, highlighting synchronization regimes and scenarios of stochastic disruption. In Section 4, we analyze intermittent chaotic behavior in detail. Finally, Section 5 summarizes our findings and discusses potential directions for future research.

2. Deterministic Model

Consider a system of two non-identical neurons modeled by the Rulkov map [16] and mutually coupled by electrical synapses:

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

(1)

Here, $f(\gamma ,u)={\displaystyle \frac{\alpha }{1+u^2}}+\gamma ,\alpha =4.1.$ In system (1), ${\gamma }_2={\gamma }_1+\Delta ,$ where $\Delta$ is a mismatch parameter. The parameter $\sigma$ regulates the strength of coupling.

2.1. Dynamics of isolated neuron

Figure 1 shows the $\gamma$-interval in which the isolated neuron $x_{t+1}={\displaystyle \frac{\alpha }{1+x_t^2}}+\gamma$ demonstrates the transformation of a stable 3-cycle into a chaotic attractor when $\gamma$ passes the crisis bifurcation point ${\gamma }_{cr}=-1.74722.$ This transition from order to chaos is justified by the largest Lyapunov exponent $\Lambda$ in the bottom panel.

2.2. Dynamics of coupled neurons

Now, consider now the corporative dynamics of the coupled system (1) with ${\gamma }_1=-1.75$ and ${\gamma }_2=-1.748$, i.e., $\Delta =0.002$. While isolated neurons with these parameters exhibit 3-cycles, the joint dynamics is governed by the coupling parameter $\sigma$, as illustrated in Figure 2, which presents the bifurcation diagram of the x variable.

At weak coupling strengths, the system (1) exhibits tristability, with three coexisting 3-cycles: $\mathcal{B}$ (blue), $\mathcal{G}$ (green), and $\mathcal{O}$ (orange). These attractor and their basins of attraction are illustrated in Figure 3 for different values of the coupling parameter $\sigma$.

Figure 3(a) illustrates the coexistence of three 3-cycles at $\sigma =0.003$. As the coupling increases, cycle $\mathcal{B}$ is annihilated, reducing system (1) to a bistable state with two coexisting 3-cycles, $\mathcal{G}$ and $\mathcal{O}$, at $\sigma =0.01$ (Figure 3(b)). A further increase in $\sigma$ eliminates cycle $\mathcal{G}$, making the system monostable with a single remaining attractor, cycle $\mathcal{O}$, at $\sigma =0.015$ (Figure 3(c)). Finally, at the crisis bifurcation occurring at ${\sigma }_{*}=0.019011$, cycle $\mathcal{O}$ loses stability, causing the system transition into chaos. The resulting chaotic attractor is shown in Figure 3(d) for $\sigma =0.05$.

The cycles $\mathcal{B}$, $\mathcal{G}$, and $\mathcal{O}$ correspond to different synchronization patterns. Solutions originating from the basins of $\mathcal{B}$ and $\mathcal{G}$ exhibit synchronization in a one-step shift mode. In the basin of $\mathcal{B}$, solutions stabilize to a 3-cycle where the y-coordinate is one step ahead of the x-coordinate, which we denote as the $[-]$ synchronization mode. Conversely, in the basin of $\mathcal{G}$, solutions stabilize to a 3-cycle for where the y-coordinate lags one step behind the x-coordinate, referred to as the $[+]$ synchronization mode. Finally, solutions in the basin of $\mathcal{O}$ stabilize to a 3-cycle without any coordinate shift, denoted as the $\left[0\right]$ synchronization mode. Thus, as the coupling increases, the system first losses the $[-]$ synchronization mode, followed by the $[+]$ synchronization mode.

3. Stochastic Model

Let us now study how noise influences system (1). The dynamics of the coupled neurons under random fluctuations in the coupling parameter, $\sigma :\sigma \to \sigma +\epsilon {\xi }_t$, is governed by the following equations:

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

(2)

Here, ${\xi }_t$ represents uncorrelated white Gaussian noise with intensity $\epsilon$.

Next, we analyze how the solutions of system (2) depend on the noise intensity $\epsilon$.

3.1. Stochastic Deformations in the Tristability Case

First, consider system (2) with $\sigma =0.003$. As shown in the previous section, at this coupling strength, the deterministic system (1) exhibits the coexistence of three 3-cycles ($\mathcal{O}$, $\mathcal{B}$, and $\mathcal{G}$) , as illustrated in Figure 3(a).

Figure 4 demonstrates how stochastic perturbations in the coupling parameter influence the system dynamics. In Figure 4(a,b,c), we plot the evolution of the difference between the x- and y-coordinates ($x-y$) of solutions, corresponding to cycles $\mathcal{O}$ (orange), $\mathcal{G}$ (green), and $\mathcal{B}$ (blue) under three different noise intensity values.

For weak noise ($\epsilon =0.0004$), stochastic solutions of system (2) that start at cycles $\mathcal{O}$ or $\mathcal{G}$ remain in these cycles, as shown in upper two panels of Figure 4(a). However, trajectories that begin at cycle $\mathcal{B}$ eventually switch to either cycle $\mathcal{G}$ or $\mathcal{O}$ after a transient period, as illustrated in lower two panels. Thus, weak noise disrupt cycle $\mathcal{B}$ while leaving the other cycles unaffected.

As shown in Figure 4(b), stronger noise ($\epsilon =0.004$) not only eliminates cycle $\mathcal{B}$ (lowest panel) but also destroys cycle $\mathcal{G}$ (middle panel) while cycle $\mathcal{O}$ remains unaffected (upper panel). Interestingly, solutions that start at $\mathcal{G}$ (middle panel) exhibit mutual transitions between $\mathcal{G}$ and $\mathcal{B}$ ($\mathcal{G}\leftrightarrow \mathcal{B}$), eventually stabilizing near cycle $\mathcal{O}$. Meanwhile, solutions originating from $\mathcal{B}$ undergo sequential transitions $\mathcal{B}\to \mathcal{G}\to \mathcal{O}$. Thus, regardless of the initial conditions, all transitions at this noise level ultimately lead to cycle $\mathcal{O}$, making the system monostable.

When the noise intensity increases to $\epsilon =0.1$, all three cycles lose stability. At this noise level, each cycle is disrupted by chaotic bursts, as shown in Figure 4(c). The system undergoes the following stochastic transitions: $\mathcal{O}\leftrightarrow \mathcal{CH}$, $\mathcal{G}\to \mathcal{O}\leftrightarrow \mathcal{CH}$, and $\mathcal{B}\to \mathcal{O}\leftrightarrow \mathcal{CH}$. Thus, regardless of the initial state, the stochastic system (2) transitions into an intermittent regime characterized by alternating small-amplitude stochastic oscillations near $\mathcal{O}$ and large-amplitude chaotic oscillations $\mathcal{CH}$. Random states of this chaotic regime are depicted in Figure 4(d).

Figure 5 illustrate the stochastic transformations the ($x-y$) variable as the noise intensity $\epsilon$ gradually increases, starting from deterministic cycles $\mathcal{B}$ (Figure 5(a)), $\mathcal{G}$ (Figure 5(b)), and $\mathcal{O}$ (Figure 5(c)). These cycles are shown in the background in light colors. Colored dots represent the random states of system (2) for solutions $\mathcal{B}$, $\mathcal{G}$, and $\mathcal{O}$, plotted after 1000 transient steps. The color of these dots matches the corresponding starting cycle but appears in a brighter shade.

For weak noise, the random solutions remain close to their respective starting cycles. However, as the noise intensity increases, stochastic switches between coexisting states begin to appear, as shown in the time series in Figure 4(a,b). These transitions are revealed in Figure 5 as dots between horizontal lines.

In the case of $\mathcal{B}$ (Figure 5(a)), these transitions disrupt the $[-]$ synchronization mode. As the noise intensity increases further, random states begin clustering around cycle $\mathcal{O}$ (see the time series in Figure 4(b)). This indicates that both the $[-]$ and $[+]$ synchronization modes are destroyed, leaving only the $\left[0\right]$ synchronization mode intact. When the noise becomes sufficiently strong, solutions starting at $\mathcal{B}$ transition directly into the intermittent regime $\mathcal{O}\leftrightarrow \mathcal{CH}$, characterized by alternating small-amplitude stochastic oscillations near $\mathcal{O}$ and large-amplitude chaotic oscillations $\mathcal{CH}$ (see the time series and random states in Figs. Figure 4(c,d)). Consequently, the $\left[0\right]$ synchronization mode is temporarily disrupted by random disturbances.

Stochastic transformations of solutions starting from the 3-cycle $\mathcal{G}$ are shown in Figure 5(b). For noise intensities up to $\epsilon ={10}^{-3}$, these solutions remain localized near the deterministic cycle. At $\epsilon =0.004$, stochastic transformations following the sequence $\{\mathcal{G}\leftrightarrow \mathcal{B}\}\to \mathcal{O}$ are observed, as illustrated in the time series in Figure 4(b) (green trace). For stronger noise, the stochastic system transits to the regime $\mathcal{G}\to \{\mathcal{O}\leftrightarrow \mathcal{CH}\}$ (see the time series shown in green in Figure 4(c) for $\epsilon =0.1$). For solutions starting at the 3-cycle $\mathcal{O}$, only a single-stage stochastic transformation occurs, where noisy oscillations near $\mathcal{O}$ transition into the intermittent mode $\mathcal{O}\leftrightarrow \mathcal{CH}$ (see the orange time series in Figure 4(c) for $\epsilon =0.1$).

To summarize the effect of noise in the tristability case, increasing noise intensity first destroys cycle $\mathcal{B}$, followed by the destruction of cycle $\mathcal{G}$, leaving cycle $\mathcal{O}$ as the most resilient to noise.

3.2. Confidence Domain Method

The cycle breakdown process described above can be analyzed analytically using the confidence domain method [26], where the mutual arrangement of confidence domains and attractors’ basins of attraction plays a key role. Below, we briefly present details of the confidence domain method.

Consider a deterministic system with an exponentially stable 3-cycle consisting of successive states $({\overline{x}}_1,{\overline{y}}_1)$, $({\overline{x}}_2,{\overline{y}}_2)$, and $({\overline{x}}_3,{\overline{y}}_3)$. Under random disturbances, the stochastic sensitivity of these states is determined by the corresponding matrices $W_1$, $W_2$, and $W_3$. These matrices satisfy the following algebraic equations [26]:

$$W_2=F_1{W}_1{F}_1^{\top }+Q_1,W_3=F_2{W}_2{F}_2^{\top }+Q_2,W_1=F_3{W}_3{F}_3^{\top }+Q_3.$$

Here, $F_t$ is the Jacobi matrix at the point $({\overline{x}}_t,{\overline{y}}_t)$ of the 3-cycle, while the matrices $Q_t$ characterize the stochastic disturbances.

For system (2), the matrix $Q_t$ is given by

$$Q_t={\left({\overline{y}}_t-{\overline{x}}_t\right)}^2\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right].$$

The stochastic sensitivity matrix $W_1$ satisfies the matrix equation

$$W_1=F{W}_1{F}^{\top }+Q,$$

where $F=F_3{F}_2{F}_1,Q=Q_3+F_3{Q}_2{F}_3^{\top }+F_3{F}_2{Q}_1{F}_2^{\top }{F}_3^{\top }.$ The other stochastic sensitivity matrices, $W_2$ and $W_3$, can be determined using the recurrence relations provided earlier.

These stochastic sensitivity matrices are useful for approximating the dispersion of random states around the points of the deterministic 3-cycle. This approximation is represented in the form of confidence ellipses. For each state $({\overline{x}}_t,{\overline{y}}_t)$ of the 3-cycle, the confidence ellipses are defined as

$${\displaystyle \frac{{\alpha }_1^2}{{\mu }_{1,t}}}+{\displaystyle \frac{{\alpha }_2^2}{{\mu }_{2,t}}}=-2{\epsilon }^{2}ln(1-P),(t=1,2,\dots ,k).$$

Here, ${\alpha }_1,{\alpha }_2$ are the coordinates of the ellipse in the basis of normalized eigenvectors $u_{1,t},u_{2,t}$ of the matrix $W_t$, while ${\mu }_{1,t},{\mu }_{2,t}$ are the corresponding eigenvalues. The parameter P represents the fiducial probability.

In Figure 6, we present magnified fragments of the phase plane of the deterministic system (1), highlighting the basins of attraction near specific states of the 3-cycles. The basin points are color-coded: light blue for cycle $\mathcal{B}$, light green for cycle $\mathcal{G}$, and white for cycle $\mathcal{O}$. A notable feature of these basins is the distinction between solid and fractal parts. The solid region surrounds the cycle points, while outside tis area, the basins of all three 3-cycles are fractally interwoven.

In addition to these deterministic basins, Figure 6 also displays confidence ellipses for the stochastic system (2). As the noise intensity increases, these ellipses expand in size. As long as the a confidence ellipse remains within the solid region of a basin, no transitions occur. However, once the ellipse extends beyond the solid part and overlaps with the fractal region, it begins to encompass points from the basins of other coexisting attractors. This expansion marks the breakdown of the initial oscillatory regime and signals transitions to alternative oscillatory modes.

Now, let us examine how the confidence domain method applies in the case of tristability.

In Figure 6(a), confidence ellipses are plotted around a point of cycle $\mathcal{B}$ for two noise intensities: $\epsilon =0.0001$ and $\epsilon =0.0004$. For $\epsilon =0.0001$, the ellipse remains entirely within the solid region of the basin, preventing any noise-induced transitions. However, at $\epsilon =0.0004$, the ellipse extends into the fractal region of the basin, leading to the stochastic destruction of oscillations near $\mathcal{B}$.

A different scenario is observed near cycle $\mathcal{G}$, where noise with $\epsilon =0.0004$ is unsufficient to disrupt stochastic oscillations. As shown in Figure 6(b), the confidence ellipse remains fully contained within the solid part of the $\mathcal{G}$ basin, ensuring stability. The same holds for cycle $\mathcal{O}$, as depicted in Figure 6(c). Consequently, a higher noise intensity is required to destabilize cycle $\mathcal{G}$. This effect is evident in Figure 6(b), where the confidence ellipse for $\epsilon =0.004$ extends into the fractal region, signaling to onset of stochastic destruction.

For cycle $\mathcal{O}$, it remains stable even under stronger noise with $\epsilon =0.05$ (see blue confidence ellipse in Figure 6(d)). However, as shown in Figure 6(d), when the noise intensity increases to $\epsilon =0.3$, the confidence ellipse extends into the fractal part of the basin, indicating the breakdown of noisy oscillations near $\mathcal{O}$.

Notably, the confidence ellipse method predicts that solutions starting from $\mathcal{B}$ and $\mathcal{G}$ transition to $\mathcal{O}$ at $\epsilon =0.004$ (see Figure 6(a,b)). In contrast, the solutions originating from $\mathcal{O}$ remain near $\mathcal{O}$ even up to $\epsilon =0.05$ (see Figure 6(d)). This reveals a broad range of noise intensity, or an "$\epsilon$-window", where all solutions of the stochastic system (2) reside near $\mathcal{O}$. This finding is supported by the results of numerical simulations presented in Figure 5. Within this $\epsilon$-window, regardless of the initial conditions, the $\left[0\right]$ synchronization mode persists.

It is worth noting that the predictions made using the confidence ellipse method align well with the results of direct numerical simulations discussed in Section 3.1.

3.3. Stochastic Deformations in the Bistability Case

We now consider system (2) with $\sigma =0.01$, which exhibits bistability in the absence of noise, characterized by the coexistence of 3-cycles $\mathcal{O}$ and $\mathcal{G}$ (see Figure 3(b)). Possible scenarios of noise-induced transformations in the distribution of random state are shown in Figs. Figure 7(a,b). In Figure 7(a), random states of solutions originating from cycle $\mathcal{G}$ are shown as a function of noise intensity $\epsilon$. A clear transition sequence is observed $\mathcal{G}\to \mathcal{O}\to \mathcal{CH}$, indicating the progressive destabilization of the initial cycle as noise increases.

The mutual arrangement of basins and confidence ellipses for $\epsilon =0.0002$ (blue) and $\epsilon =0.001$ (red) is illustrated in Figure 7(c). For $\epsilon =0.0002$, the confidence ellipse lies entirely within the solid part of the basin of $\mathcal{G}$, meaning that random solutions exhibit only slight fluctuations around $\mathcal{G}$. However, for $\epsilon =0.001$, the ellipse extends into the fractal region of the basin, allowing for noise-induced transitions $\mathcal{G}\to \mathcal{O}$. The subsequent transition $\mathcal{O}\to \mathcal{CH}$, observed in Figs. Figure 7(a,b), is accurately predicted by the confidence ellipses shown in Figure 7(d).

3.4. Stochastic Deformations in the Monostability Case

Let us now examine the behavior of the stochastic system (2) in the case of monostability, when the original deterministic system (1) possesses only one attractor: the in-phase 3-cycle $\mathcal{O}$. This cycle, corresponding to $\sigma =0.015$, is shown in Figure 3(c). For this parameter value, the evolution of random states originating from $\mathcal{O}$ is presented in Figure 8(a) as a function of noise intensity $\epsilon$.

As evident from the figure, the 3-cycle $\mathcal{O}$ demonstrates notable resilience to noise. However, beyond a certain of threshold intensity of $\epsilon$, the system (2) undergoes an abrupt change in the probability distribution of its states. Specifically, in Figure 8(b), the random states for $\epsilon =0.02$ (blue) remain tightly clustered near the deterministic cycle $\mathcal{O}$, whereas for $\epsilon =0.1$ (orange), the states are broadly scattered across the phase space. This behavior reflects a clear manifestation of noise-induced excitation.

While in the cases of multistability, sharp deformations in the distributions were attributed to noise-induced transitions between coexisting attractors, the mechanism driving the transformation in this monostable case is different. Here, the influence of transient attractors in the underlying deterministic system (1) becomes significant.

4. Transients and Intermittent Synchronization

In this section, we investigate the transient dynamics of the deterministic system (1) and explore noise-induced intermittent synchronization in the corresponding stochastic system (2).

4.1. Short and Long Transients Basins

Although the monostable system (1) with $\sigma =0.015$ ensures that all trajectories eventually converge to the 3-cycle $\mathcal{O}$, the rate of convergence strongly depends on the initial conditions. This behavior is illustrated in Figure 9(a), where initial points of solutions that take more than 200 iterations to converge to $\mathcal{O}$ with an accuracy of 0.001 are marked in orange. These points define the long transients basin (LTB), while the remaining region—corresponding to faster convergence—is marked in white and referred to as the short transients basin (STB). Notably, the boundary between these basins exhibits a complex, fractal structure. In Figure 9(b), an enlarged fragment of the phase plane is shown near one of the points of 3-cycle $\mathcal{O}$, highlighting both the STB and LTB. The STB consists of a solid (white) region and a fractal part, represented by a mixture of white and orange areas. The selected 3-cycle point lies within the solid part of the STB.

We now examine how the transient behavior of solutions to the deterministic system (1) depends on the choice of the initial condition. Specifically, we consider two nearby points, $A(2.32,2.28)$ and $B(2.32,2.27)$. Point A lies within the solid part of the STB, whereas point B belongs to its fractal region.

During transients, the system may alternate between synchronous ($x=y$) and asynchronous ($x\ne y$) states. To quantify these states, we introduce the synchronization index defined as

$$z_t=sgn\left[(x_{t+1}-x_t)(y_{t+1}-y_t)\right].$$

Note that $z_t=1$ when the neurons evolve synchronously (i.e., both variables increase or decrease together), and $z_t=-1$ otherwise. Examples of synchronous and asynchronous transients are shown in Figs.Figure 9(c) and Figure 9(d), respectively. Specifically, Figure 9(c) illustrates that the trajectory starting at point A converges to the 3-cycle $\mathcal{O}$ while maintaining in-phase synchronization ($z_t=1$). In contrast, Figure 9(d) shows that the trajectory from point B exhibits a large-amplitude burst with temporary desynchronization, characterized by the interval where $z_t=-1$.

The diversity of transient behaviors for various initial conditions $(x_0,y_0)$ near point B—located in the fractal region—is illustrated in Figure 9(e). Notably, for $y_0=2.2712$, the system exhibits three bursts before ultimately converging to the 3-cycle $\mathcal{O}$.

The underlying reason for the noise-induced excitability described above and illustrated in Fig.Figure 8 lies in the presence of a fractal region associated with large-amplitude burst transient attractors in the deterministic system (1). For weak noise, random trajectories originating near the cycle $\mathcal{O}$ remain close to its states within the solid part of the short transient basin (STB), resulting in only minor stochastic fluctuations around the cycle. As the noise intensity increases, these trajectories may cross into the fractal region, triggering large-amplitude burst oscillations. This qualitative shift is anticipated by the positions of the confidence ellipses. In Figure 9(b), the ellipse for $\epsilon =0.02$ lies entirely within the solid part of the STB. However, the ellipse for $\epsilon =0.1$ partially overlaps the "dangerous" fractal zone. This overlap serves as an indicator of the onset of noise-induced excitation. The prediction aligns well with the results of direct numerical simulations (compare Figs. Figure 9(b) and Figure 8).

It is important to note that the distinction between the solid and fractal regions in Figs.Figure 9(a,b) is based on a fixed threshold of $N=200$ iterations. However, the convergence rate of solutions of system (1) to cycle $\mathcal{O}$ depends on the specific initial condition. This dependence is illustrated in Figure 10, where the color map represents the number of iterations N required to reach the cycle with a convergence accuracy of 0.001. The contrasting characteristics of the solid and fractal regions are clearly visible in this representation. Notably, the confidence ellipse for $\epsilon =0.1$ intersects the fractal region, indicating the potential for noise-induced excitation.

4.2. Intermittent Synchronization

We now examine how the synchronous behavior of the stochastic system (2) evolves in the excitability regime as the noise intensity increases. Figure 11 presents time series of the variable x and the synchronization indicator z for solutions of the stochastic system (2), initialized at the deterministic 3-cycle $\mathcal{O}$. Two noise intensities are considered: $\epsilon =0.05$ and $\epsilon =0.1$.

The time series reveal two distinct phases: laminar and turbulent. During the laminar phase, the system maintains in-phase synchronization ($z=1$), with trajectories fluctuating slightly around the points of the periodic 3-cycle. In contrast, the turbulent phase is characterized by burst-like oscillations and a breakdown of synchronization. It is important to note that the durations $\tau$ of the laminar intervals are random and tend to decrease as the noise intensity increases—compare Figure 11(a) and Figure 11(b).

Since the turbulent phase emerges due to the presence of noise, this phenomenon belongs to noise-induced intermittency. Similar behavior has been observed in previous studies on systems with randomly modulated parameters [29,30,31,32,33]. Intermittency is typically characterized by the distribution of laminar phase durations, computed for fixed control parameters, and by analyzing the dependence of the mean laminar duration on a critical parameter.

Figure 12(a) illustrates the dependence of the mean duration of laminar phase $\langle \tau \rangle$ on noise intensity $\epsilon$ for different values of the coupling strength $\sigma$. For $\sigma =0.012$, $\sigma =0.015$, and $\sigma =0.017$, the deterministic system (1) possesses a single in-phase 3-cycle attractor, whereas for $\sigma =0.02$, the system exhibits a single chaotic attractor. The sharp decline in $\langle \tau \rangle$ for the 3-cycle cases marks the onset of noise-induced excitation, characterized by the emergence of temporal bursts. In contrast, the horizontal plateau observed in the chaotic regime suggests that the mean interburst intervals remain unchanged as noise intensity increases.

Figure 12(b) shows the probability of the stochastic system (2) being in the co-directed mode as a function of noise intensity $\epsilon$ for various coupling strengths $\sigma$. For all considered values of $\sigma$, including the chaotic regime ($\sigma =0.02$), the system predominantly remains in the co-directed mode, with probabilities $P>0.5$.

Figure 12(c) presents the dependence of the mean laminar duration $\langle \tau \rangle$ on the deviation parameter d for different values of $\epsilon$. As expected, the farther $\sigma$ is from the crisis bifurcation point ${\sigma }_{*}$, the longer the mean duration of laminar (co-directed) behavior.

Finally, in Figure 13, we present the largest Lyapunov exponent $\Lambda$, which serves as a key indicator of system stability. This exponent provides a quantitative criterion for distinguishing between regular dynamics ($\Lambda <0$) and chaotic behavior ($\Lambda >0$).

Figure 13(a) shows $\Lambda$ as a function of noise intensity $\epsilon$ for various coupling strengths $\sigma$, while Figure 13(b) illustrates the dependence of $\Lambda$ on both control parameters, $\sigma$ and $\epsilon$. From these graphs, one can estimate the threshold noise intensity at which the system transitions from order to chaos. For instance, at $\sigma =0.015$, this transition occurs around $\epsilon \approx 0.03$. It is worth noting that this noise level also corresponds to the onset of stochastic excitation, as previously illustrated in Figure 8(a).

5. Conclusion

In this study, we examined the stochastic disruption of synchronization in a system of two non-identical Rulkov neurons coupled via an electrical synapse. Our analysis uncovered regions of mono-, bi-, and tristability, demonstrating how distinct synchronization regimes arise as a function of coupling strength. Through parametric exploration, we identified conditions under which stochastic perturbations to the coupling parameter induce a transition from lag synchronization to monostable in-phase synchronization.

A central finding of this work is the noise-induced intermittency within the synchronous 3-cycle regime, where increasing noise intensity leads to transient episodes of desynchronization characterized by chaotic bursts. We showed that this behavior is closely linked to the presence of chaotic transient basins, which act as precursors to the disruption of stable synchronization. Moreover, we identified a specific range of noise intensities where this intermittency persists, with its statistical properties—such as laminar phase duration—being highly sensitive to the coupling strength.

Methodologically, the study combines direct numerical simulations with an analytical framework based on confidence ellipses, providing a powerful tool for identifying and predicting the onset of noise-induced excitations. These findings advance our understanding of how stochastic effects influence synchronization in coupled neuronal systems and offer valuable insights into noise-driven dynamics in both biological and artificial neural networks.

Future directions include extending this framework to larger and more complex networks, exploring the role of network topology, and investigating practical implications for neuromorphic engineering and biomedical modeling of brain activity under noisy conditions.