Delay-induced multirhythmicity in a 3-cell inhibitory network

Introduction

Central pattern generators (CPGs) are neural circuits that produce rhythmic motor patterns without rhythmic external input [1, 2, 3]. Interactions within biological CPGs are not instantaneous: finite synaptic transmission and signal-propagation times introduce delays into network coupling [4, 5, 6]. Such delays can substantially reorganize the network dynamics, affecting synchronization, phase locking, and the coexistence of rhythmic states even in small networks of coupled relaxation oscillators [7, 8, 9]. To analyze multistability in coupled CPG networks, Shilnikov and collaborators developed a return-map framework that provides a geometric description of phase-locked rhythms and transitions among them in small neuronal motifs [10]. This approach has proved particularly effective for identifying coexisting rhythms, determining their basins of attraction, and characterizing the global organization of the network dynamics [11, 12, 13, 14, 15]. In our previous work [7], we extended this framework to delayed inhibitory networks, enabling a systematic analysis of how synaptic delays affect rhythm generation, multistability, and transitions among rhythmic states. Phase reduction provides a complementary description of coupled-oscillator dynamics. In the weak-coupling regime, a high-dimensional network of oscillators can be reduced to a low-dimensional system governing their phase interactions [16, 17, 18, 19]. Early studies of delay-coupled oscillators incorporated transmission delays explicitly into the phase arguments [20, 21, 22]. Subsequent work established that, under weak coupling, these delays can instead be represented as effective phase shifts in the interaction functions of the reduced equations [9, 16, 23, 24, 25, 26, 27] as long as the delay is not too long [28].

Previous studies have shown that, in CPG networks without delays, phase-difference models accurately reproduce the fixed-point organization and stability properties revealed by return maps [19]. Whether this correspondence persists in delayed CPG networks, however, remains largely unexplored. Delays can generate more complex dynamical structures, including invariant curves, whose location and stability are not readily determined from return-map constructions alone. This limitation motivates a systematic comparison of return-map and phase-reduced descriptions of delayed CPG networks.

In this work, we investigate a delayed three-cell inhibitory CPG network by combining return-map analysis with phase reduction of a network of generalized FitzHugh–Nagumo oscillators [7, 12, 14, 15]. Using phase-response curves (PRCs) and coupling functions [16], we derive a two-dimensional phase-difference system and systematically compare its predictions with the dynamics of the full delayed network. We show that the phase-reduced model not only reproduces the principal rhythmic structures revealed by the return maps but also provides a more detailed characterization of the phase space. In particular, it identifies stable synchronous states with extremely small basins of attraction that can be difficult to detect through direct return-map simulations. The reduced model also permits the location and stability of limit cycles in the phase-difference plane to be determined directly, thereby clarifying the dynamical origin of the invariant curves observed in the return maps. Using both the full and reduced models, we perform extensive two-parameter analyses of rhythm generation and transitions among coexisting states. The close agreement between the two descriptions establishes phase reduction as an accurate and computationally efficient framework for large-scale exploration of delayed CPG networks.

Methods

We consider three weakly coupled oscillators arranged in a symmetric networking ring, with a transmission delay in each coupling term. The state of the $j$-th oscillator, denoted by $\mathbf{X}_j(t) \in\mathbb{R}^n$, evolves according to

where

Here, $F:\mathbb{R}^n \to\mathbb{R}^n$ describes the intrinsic dynamics of each oscillator, $\tau\ge0$ denotes the transmission delay, and $0 < g_{j-1,j} \ll1$ are small coupling parameters. The functions $G:\mathbb{R}^n \times\mathbb{R}^n \to\mathbb{R}^n$ specifies the interaction between a postsynaptic oscillator and a delayed presynaptic oscillator. To identify activation events, we introduce a scalar observable

where $\mathcal{V}:\mathbb{R}^n \to\mathbb{R}$ maps the state of each oscillator to the scalar component or output variable used to detect threshold crossings.

Time delays are an intrinsic feature of neural signal transmission in biological systems, and their inclusion in CPG models is essential for accurately capturing signal propagation and enhancing biological plausibility. Let $T_0$ denote the period of the uncoupled oscillator. Then, we can normalize the transmission delay $\tau$ by this intrinsic period as following: $\alpha=\tau/T_0$. We will assume that $0 \le\alpha\le1$.

For a periodic phase-locked solution of the full network, let $T$ denote its common period. On a $T$-periodic network solution, $\mathbf{X}_j(t-kT)=\mathbf{X}_j(t)$ for every integer $k$. Thus, for $\tau=kT\ge0$, the delayed input coincides with the instantaneous input along that solution. This identity does not imply identical stability or transient dynamics in the delayed and undelayed networks. Noninteger-period delays generally modify the relative timing of the inputs and can reorganize the rhythmic states.

Return map

Because the full higher-dimensional network comprises three coupled oscillatory cells, its relative phase dynamics can be described by two independent phase-difference variables, $\Delta=(\Delta_{12},\Delta_{13})$. Here, $\Delta_{1j}$ ($j=2,3$) denotes the phase lag of node $j$ relative to node 1. Together, these variables constitute the maximal set of independent phase differences among the three oscillators. Following the procedure of Wojcik et al. [13], we obtain the phase lags from threshold-crossing times of the membrane potentials. Specifically, we choose node 1 as the reference oscillator and record the sequence of times $t_1^{(k)}$ at which $V_1(t)$ crosses the threshold $V_1=0$ from below, with each crossing corresponding to the $k$th activation of node 1. For each reference event, we then identify the subsequent activation times $t_j^{(k)}$, $j=2,3$, of nodes 2 and 3. The time lags of these events relative to the reference event are normalized by the corresponding reference period $T_1^{(k)}=t_1^{(k+1)}-t_1^{(k)}$. This procedure yields a sequence of phase-lag pairs along the forward trajectory, denoted by $\Delta^{(k)}=(\Delta_{12}^{(k)},\Delta_{13}^{(k)})$, where $\Delta_{12}^{(k)}=(t_2^{(k)}-t_1^{(k)})/T_1^{(k)}$, $\Delta_{13}^{(k)}=(t_3^{(k)}-t_1^{(k)})/T_1^{(k)}$. The cycle-to-cycle evolution of the phase lags is governed by the Poincaré map $\Pi:P_k\mapsto P_{k+1}$, where $P_k=\Delta^{(k)}$. This map advances the phase-lag vector from one activation cycle of the reference node to the next. Its dynamics determine whether trajectories converge to stable fixed points (FPs), corresponding to phase-locked rhythms with constant phase lags, or exhibit more complex behavior, such as motion on an invariant closed curve in the phase-lag plane. Stable FPs represent asymptotically stable phase-locked rhythms, whereas the coexistence of multiple attractors indicates network multistability.

When represented on the unit square with both coordinates defined modulo 1, the two-dimensional phase-lag map preserves the topology of a torus, with opposite edges identified. At the pairwise level, the basic phase relationships include $\Delta=0$ for cells in-phase and $\Delta=0.5$ for cells in anti-phase. To identify the attractors and their basins of attraction, we initialize the network on a grid of phase-lag values (e.g., $20 \times20$) and compute sufficiently long sequences of network activity for the transient dynamics to decay. Fixed points with constant $(\Delta_{12}^{*},\Delta_{13}^{*})$ lags correspond to stable phase-locked rhythms. Initial conditions assigned the same color converge to the same attractor and therefore belong to the same basin of attraction. Distinct rhythmic states are represented by fixed points at different locations in the phase-difference plane.

We classify the observed rhythms according to the locations and stability of the fixed points (FPs) and the geometry of the asymptotic trajectories in the phase plane. Traveling-wave (TW) rhythms are associated with stable FPs near $(\Delta_{12}, \Delta_{13}) \approx(1/3, 2/3)$ and its symmetry-related counterpart. Pacemaker (PM) rhythms correspond to stable FPs near $(0, 1/2)$, $(1/2, 1/2)$, and $(1/2, 0)$, whereas synchronous (SY) rhythms are identified by convergence to $(\Delta_{12}, \Delta_{13}) \approx(0, 0)$. Trajectories that converge to a closed invariant curve in the phase plane are associated with the phenomenon of phase jitter – slow oscillations centered around a phase-locked state, whereas trajectories exhibiting persistent phase drift are classified as phase-slipping (PS) dynamics.

Standard phase reduction

We assume that each uncoupled oscillator possesses an exponentially asymptotically stable limit cycle $\hat{X}(t)$ with period $T_0$. Let $\hat{Z}(t)$ denote the corresponding $T_0$-periodic infinitesimal phase-response curve (PRC), normalized according to

For sufficiently weak coupling, the delayed three-cell network introduced above can be reduced to a phase model [16].

The interaction function $H$ is defined by

Because the phase variables are defined modulo 1 and $\hat{X}$ is $T_0$-periodic, the interaction function $H$ is 1-periodic. The PRC determines the phase sensitivity of each oscillator to weak coupling perturbations. For an oscillator at a given phase, as illustrated by the colored circles in Figure 1(b), each delayed coupling input produces a small phase correction whose magnitude and direction depend on the phase at which the perturbation is received. Through repeated interactions, these phase corrections continuously modify the relative phases between the oscillators. Depending on the initial phase and coupling function, the network can exhibit different stable phase-locked rhythms, with representative examples shown in Figure 1(c).

Denoting the normalized phase of oscillator $i$ by $\phi_i \in [0, 1), to leading order in the weak-coupling limit, we obtain [16]

where $\Omega= 1/T_0$ and $g_{ji}$ is the coupling strength from oscillator $j$ to oscillator $i$. Since one uncoupled oscillation corresponds to one unit of phase, the transmission delay $\tau$ corresponds to the phase delay

where $T_0$ is the period of the uncoupled oscillator and $\Omega= 1/T_0$. When $\alpha$ is small compared to $1/(\max g_{ij})$ it enters the reduced equations as a fixed phase shift in the interaction function [28]. The resulting system (6) consists of three phase equations; eliminating the common phase yields a two-dimensional system governing the two independent phase differences.

Choose phase zero at the reference threshold crossing, so that $\varphi_1(t_1)=0$ modulo one. For uncoupled oscillators the time-to-phase conversion is exact with period $T_0$. Under weak coupling it holds to leading order for transverse crossings and bounded inter-event intervals:

Hence, defining the reduced phase differences by

gives the same sign convention as the event-based lags. When $T_1^{(k)}=T_0+O(g)$, normalization by the measured reference interval instead of $T_0$ also changes these lags only by $O(g)$. The reduced phases and the measured event lags are thus leading-order counterparts, rather than exactly identical coordinates. These variables evolve at rates of order $g$, and hence on a slow time scale of order $g^{-1}$. To leading order in $g$, they satisfy

The fixed points (FPs), given by the intersections of the nullclines of Eq. (9), correspond to phase-locked solutions of the original oscillator network. They are located at the intersections of the two nullclines, which are defined by the conditions: $\dot{\Delta}_{12}=0$ and $\dot{\Delta}_{13}=0$, and are shown as the red and blue curves, respectively, in Figure 2. Linearization of Eq. (9) about a FP $\Delta^*$ yields the Jacobian matrix $J$. The eigenvalues of $J$ determine the local stability of the corresponding phase-locked rhythm.

Numerics

Simulations of the full network equations were performed with a fixed time step $dt=0.01$. The instantaneous state variables were advanced using a classical fourth-order Runge–Kutta scheme, while the delayed presynaptic terms were evaluated from a discrete history buffer. For $t\in[-\tau,0)$, the delayed presynaptic variables were initialized to zero, whereas the oscillator states at $t=0$ were specified by the chosen initial conditions. Within each integration step, the delayed input obtained from the history buffer was held fixed across the four Runge–Kutta stages, without additional interpolation. Step-size convergence was checked by repeating the simulations with a smaller time step, with no appreciable changes in the reported dynamical behavior. The resulting trajectories are processed using Python scripts for visualization and return-map construction. The nullclines of the phase-difference system obtained using MATLAB were double-checked with calculations using XPPAUT [29]. Through systematic parameter variation, we track changes in the geometry of the return map and identify bifurcations and transitions among the network rhythms. This analysis reveals dynamical structures that are not readily inferred from voltage traces alone and characterizes the parameter dependence, stability, and multistability of the network-generated rhythms.

Results

Generalized FitzHugh–Nagumo network

To compare and validate the return-map and phase-reduction approaches, we use a generalized FitzHugh–Nagumo (FHN) model as a concrete realization of each neuronal oscillator. The FHN model is a canonical low-dimensional reduction of excitable membrane dynamics and can be viewed as a biologically motivated extension of the van der Pol relaxation oscillator. It retains the essential interaction between fast voltage activation and slow recovery while omitting the detailed ionic conductances of Hodgkin–Huxley-type models. The generalized FHN model used here incorporates an explicit slow–fast time-scale separation and a sigmoidal recovery nullcline, allowing it to reproduce relaxation oscillations and the envelope-level alternation between active and quiescent phases observed in neuronal bursters, both in isolation and under synaptic perturbations [12]:

Here, the state of the $i$th neuron is described by the fast membrane-voltage variable $V_i$ and some slower recovery or voltage-dependent gating $h_i$-variable. The small parameter $0 < \varepsilon< 1$ controls the separation between the fast voltage and slow recovery time scales, whereas $I$ is an applied current. The parameter $V_0$ determines the midpoint of the sigmoidal $h$-nullcline, and $k$ controls its steepness. The cubic $V$-nullcline is determined by the voltage equation and depends on $I$ and the synaptic input. Unless stated otherwise, we use $k = 10$ and $V_0 = 0$. We also set the synaptic activation threshold to $V_{\mathrm{th}} = 0$. The time-scale parameter $\varepsilon$ should not be confused with the weak-coupling parameter $g$ introduced in the general network formulation.

A presynaptic neuron activates its outgoing synapse when its voltage exceeds $V_{\mathrm{th}}$. Synaptic interactions are modeled using the fast threshold modulation (FTM) formalism,

Here, $\tau$ is the transmission delay and is treated as a bifurcation parameter. The steep sigmoidal function $\Gamma$ approximates a binary synaptic switch: it is close to 1 when the delayed presynaptic voltage exceeds $V_{\mathrm{th}}$ and close to 0 when it lies below the threshold. Thus, the FTM formalism distinguishes between the active (“on”) and inactive (“off”) states of the synaptic gate.

The nonnegative weight $g_{ji}$ specifies the strength of the directed synapse from neuron $j$ to neuron $i$. The sign of the synaptic current is determined by the driving force $V_{\mathrm{rev}} - V_i$. A reversal potential below the relevant voltage range, such as $V_{\mathrm{rev}} = -1.5$, produces inhibitory coupling, whereas a reversal potential above that range, such as $V_{\mathrm{rev}} = 1.5$, produces excitatory coupling. In the present study, we focus on inhibitory interactions.

It is well established [7] that delayed three-node CPG networks can generate three canonical classes of rhythms: a synchronous state (SY; yellow point), three pacemaker states (PM; green, red, and blue points), and two traveling-wave states (TW; purple and black points), as shown in Figure 2 and Figure 4. Cyan lines and cyan circles represent phase-slipping (PS) and limit-cycle (LC) dynamics, respectively, with the corresponding phase trajectories shown in gray. Purple diamonds denote saddle FPs, whereas black squares denote unstable FPs.

We compare the phase-space structures obtained from the return map and the phase model in the $(\Delta_{12},\Delta_{13})$–phase-lag map. The two approaches show good agreement across several parameter sets, and one representative example is shown in Figure 2. The locations and stability of the FPs inferred from the return map agree closely with those determined from the nullcline intersections for locating FPs, while their stability is determined by local linearization of the phase model. Direct integration of the phase-difference system also identifies stable LCs corresponding to invariant closed curves. In some parameter regimes, the basin of attraction of the synchronous state at the origin is extremely small, making

the stable yellow SY FP difficult to detect from return-map trajectories alone. The phase model resolves this difficulty by determining the stability of the corresponding FP directly from the eigenvalues of the linearized phase-difference system, independently of its basin size.

Using both return-map analysis and phase reduction, we systematically investigate the dynamics of the delayed three-neuron network across different parameter regimes. The return-map framework is used to analyze the effects of coupling strength and transmission delay directly in the full network. The corresponding phase model provides the precise locations and stability of saddle and unstable FPs and identifies the LCs associated with invariant curves in the return map. We also use the reduced model to explore the $(\varepsilon,I)$ parameter plane, for which its lower dimensionality permits denser parameter sweeps and a more detailed classification of the resulting rhythms at substantially lower computational cost. Together, the two approaches reveal a rich organization of coexisting rhythmic attractors and the bifurcations that mediate transitions among them.

Return-map analysis of asymmetric inhibitory networks

We investigate how asymmetry in the inhibitory coupling affects the rhythmic dynamics of the delayed three-neuron network. Specifically, the reciprocal coupling pair $g_{13} = g_{31} \equiv g_{13/31}$ is varied relative to the remaining synaptic weights. Using Poincaré return maps together with a systematic parameter sweep, we classify the multirhythmic dynamics in the two-parameter, $(g_{13/31},\alpha)$, plane.

Figure 3 shows the global dynamical organization of the delayed three-neuron motif in the $(g_{13/31},\alpha)$-parameter plane. The distribution of dynamical regimes is approximately reflected about $\alpha=0.5$, although this symmetry is not exact. Because $\alpha$ and $1-\alpha$ correspond to phase shifts of opposite orientation modulo one period, the approximate reflection is consistent with the exchange symmetries of the ring architecture. Its imperfect character reflects the detailed shape of the interaction function and the nonlinear dynamics of the full network.

Pacemaker rhythms dominate for both short and nearly full-period delays, particularly for $\alpha\in(0,0.1)$ and $\alpha\in(0.9,1)$. In contrast, PS dynamics are distributed

approximately symmetrically about $\alpha=0.5$. Representative parameter regions and their characteristic rhythms are examined below. Unless stated otherwise, the intrinsic parameters are fixed at $I=0.6$ and $\varepsilon=0.3$ for all three neurons.

Figure 4 illustrates the Poincaré return maps for a short synaptic delay, $\alpha=0.08$, as the reciprocal coupling strength $g_{13/31}$ is varied. For weak coupling, see Figure 4(a), a stable green PM coexists with PS dynamics. Unlike a phase-locked attractor, the PS state does not converge to a single FP or invariant closed curve. Instead, trajectories initialized from different points follow distinct, non-convergent drifting paths on the phase torus, indicating persistent phase drift without synchronization or stable phase locking.

As $g_{13/31}$ increases, two stable traveling wave FPs or TWs located symmetrically near $(1/3,2/3)$ and $(2/3,1/3)$ in the map and shown in black and purple are observed in Figure 4(b). Although stable, these TWs initially possess small basins of attraction and are therefore strongly sensitive to the initial phase configuration. With a further increase in coupling, the PM and TW FPs disappear, and the dynamics become dominated by two stable LCs surrounding the former TW regions [Figure 4(c)].

When $g_{13/31}=0.001$, the coupling symmetry is restored. A stable SY state is then observed with a very small basin of attraction, leaving the phase space predominantly occupied by the two LCs [Fig. 4(d)].

Figure 5 shows the evolution of the Poincaré return maps for the asymmetric motif with $g_{13/31}=0.0006$ as the synaptic delay increases. For relatively short delays, the network dynamics are dominated by PS [Fig. 5(a)]. As the delay increases, a stable SY state is observed in Figure 5(b) and subsequently coexists with PS dynamics in Figure 5(c). Finally, when the delay approaches one full oscillation period, a green PM emerges and becomes the dominant rhythm [Fig. 5(d)].

Figure 6 illustrates multistability in the symmetric motif with $g_{ij}=0.001$. For short delays, three PM states coexist with two stable LCs [Fig. 6(a)]. In the symmetric phase model, the SY equilibrium exists for every $\alpha$-value because all three phase velocities coincide when the phases coincide. Thus an observed onset of stable synchrony must be distinguished from creation of the equilibrium. As the delay increases, stable SY and TW states are observed [Fig. 6(b)]. For longer delays, the TWs dominate the phase space [Fig. 6(c)]. As the delay approaches one oscillation period, the PMs reappear and again become dominant [Fig. 6(d)]. This sequence is qualitatively consistent with the organization reported for $g_{\mathrm{all}}=0.0005$ in Ref. [7], indicating that the principal polyrhythmic structure persists as the coupling strength is increased.

Weak-coupling analysis via phase reduction

Although the return-map approach provides a detailed description of the full network dynamics, its computational cost increases substantially in the weak-coupling regime. As the coupling strength decreases, trajectories approach their asymptotic states on a slow time scale of order $g^{-1}$. Accurate classification of the attractors and their basins therefore requires long integration times at every point in a parameter sweep, making direct return-map construction inefficient for large-scale analyses.

To address this limitation, we employ standard phase-reduction theory. Under sufficiently weak coupling, the full network is approximated by a low-dimensional system of coupled phase oscillators. To first order, the reduced equations are determined by the infinitesimal phase-response curve (PRC) and the associated interaction function. The PRC quantifies the sensitivity of an oscillator’s phase to a weak perturbation and can be determined numerically or experimentally. Once the limit cycle, PRC, and interaction function of the individual oscillator have been computed, the phase-reduction framework can be applied systematically to weakly coupled networks with different coupling architectures. The resulting phase-difference equations permit the direct calculation of phase-locked states, their stability, and the bifurcations that organize transitions among network rhythms. We apply this framework in the weak-coupling regime $g = 0.0005$ to investigate how the intrinsic time-scale parameter $\varepsilon$ and the applied current $I$ jointly affect the rhythmic dynamics. These quantities are treated as bifurcation parameters in a two-parameter analysis of the phase model.

Because changing either parameter modifies the intrinsic limit cycle, the corresponding PRC and interaction function are evaluated for each parameter pair, shown in Figure 7. The lower dimensionality of the reduced system then permits a detailed classification of stable FPs and other phase-space structures across the $(\varepsilon, I)$ parameter plane.

The most pronounced changes in rhythmic organization occur for relatively short delays. We therefore focus on three representative cases occurring specifically at $\alpha= 0.05$, $0.10$, and $0.20$, as shown in Figure 8. The resulting diagrams reveal the global organization of the SY, PM, and TW states and their coexistence in the reduced phase space. Consistent with Ref. [19], the rhythmic regimes exhibit an approximate reflection between the low- and high-current regions about an intermediate range of $I$. This organization is not exactly symmetric,

reflecting the dependence of the limit cycle and interaction function on the applied current.

At $\alpha= 0.05$, we examine the rhythmic dynamics along two parameter paths: varying $\varepsilon$ at fixed $I = 0.425$ in Figure 9, and varying $I$ at fixed $\varepsilon= 0.15$ in Figure 10. In Figure 9(a), three PMs and two TWs coexist with the colored trajectories showing their corresponding basins of attraction. As $\varepsilon$ increases, the two TWs change from stable in Figure 9(b) to unstable foci and two invariant closed curves emerge in Figure 9(c) through an Andronov–Hopf bifurcation. In the continuous phase-difference system, these invariant curves correspond to stable LCs. Specifically, the corresponding eigenvalues change from $-0.0003 \pm0.0094i$ to $0.0001 \pm0.0119i$, indicating that the real part of a complex-conjugate pair crosses the imaginary axis at approximately $\varepsilon_H = 0.255$. The same Jacobian and eigenvalue analysis was also performed for the TWs in Figure 10, Figure 11 and Figure 12, confirming the occurrence of the corresponding Andronov–Hopf bifurcations. With a further increase in $\varepsilon$, the two saddles near each PM gradually approach and merge with the right-hand saddle in a subcritical pitchfork bifurcation, causing the PM to lose stability. Consequently, the dynamics are dominated by stable LCs, while a stable SY state coexists with a small basin of attraction in Figure 9(d).

For $\varepsilon= 0.15$, two stable TWs dominate the dynamics at $I = 0.50$ [Figure 10(a)]. Then three PMs appear and coexist with two TWs in Figure 10(b). When $I$ is increased to $0.58$, the TW FPs lose stability while two LCs emerge through Andronov–Hopf bifurcations around their locations in the return map in Figure 10(c). Finally, the LCs disappear and PM states become the dominant rhythm in Figure 10(d).

Figure 11 illustrates transitions among multiple coexisting attractors at the fixed delay $\alpha= 0.10$. At $\varepsilon= 0.15$, two stable TWs dominate most of the phase space, while a smaller set of initial conditions converges to the stable SY state [Figure 11(a)]. As $\varepsilon$ increases to $0.18$, both two stable LCs emerge, thus coexisting with the three persisting attractors: two TWs and the SY state, as documented in Figure 11(b).

With a further increase in $\varepsilon$, the small unstable LCs gradually shrink toward the stable TW equilibria and eventually disappear as the TWs lose stability through subcritical Andronov–Hopf bifurcations, while the two stable LCs persist. Meanwhile, the unstable fixed points approach and eventually collide with nearby saddles, leading to their annihilation, while larger unstable LCs emerge via homoclinic saddle-node bifurcation in Figure 11(c). For fixed $\varepsilon= 0.28$, as $I$ increases to $0.43$, two TWs undergo Andronov–Hopf bifurcations and regain stability, accompanied by the emergence of two small unstable LCs surrounding them. The larger unstable LCs and the two stable LCs persist, as shown in panel (d). Now, the phase-difference system exhibits five coexisting attractors: two stable TWs, two stable LCs, and one stable SY state. The small unstable LCs separate the basins of attraction of the stable TWs and the surrounding stable LCs, while the larger unstable LCs delimit the basin of attraction of the SY state.

Finally, we examine the dynamics at $\alpha= 0.20$. For $I = 0.50$, the SY state dominates the phase space, see Figure 12(a). When the current is increased to $I = 0.58$, both repelling TWs regain stability through Andronov–Hopf bifurcations, accompanied by the emergence of two small unstable LCs surrounding them, and now coexist with the synchronous rhythm [Figure 12(b)].

Conclusion

We investigated the dynamics and transformations of delayed three-cell inhibitory CPG networks by combining return-map analysis of the full generalized Fitzhugh–Nagumo system with a corresponding phase-difference model. The two approaches provide complementary descriptions of the network dynamics. Return maps retain a direct connection to the full system and reveal the global organization of rhythmic trajectories and their basins of attraction. Within the weak-coupling approximation, the phase-difference model provides a more computationally efficient and quantitatively precise description of the relative-phase dynamics. It permits direct localization of phase-locked states through nullcline intersections, determination of their stability from the linearized system, and identification of stable limit cycles corresponding to invariant closed curves observed in the return maps.

The phase model is particularly useful when an attractor has a very small basin of attraction. Such states can be difficult to detect through direct forward simulations because only a small set of initial conditions converges to them. In contrast, the corresponding equilibria of the phase-difference system can be located and classified independently of their basin size. This distinction allowed us to identify the synchronous FP and determine its stability even in parameter regimes where its basin is extremely small.

Using the return-map framework, we first examined the joint effects of the normalized delay $\alpha$ and the reciprocal coupling strength $g_{13}/g_{31}$, varied relative to the remaining inhibitory synapses. The resulting parameter scans revealed a rich multistable landscape containing synchronous, pacemaker, traveling-wave, phase-slipping, and limit-cycle dynamics, together with numerous coexistence regimes. Phase-slipping dynamics persist across a broad region of parameter space, while variations in coupling asymmetry and delay organize transitions among synchronous, traveling-wave, and pacemaker rhythms.

We then used the phase-reduced model to perform two-parameter scans in the intrinsic time-scale parameter $\varepsilon$ and the applied current $I$ for relatively short delays. The reduced system reproduced the global organization of the principal rhythms and revealed an approximate reflection between the low- and high-current regimes. This reflection is accompanied by an approximate reversal of the PRC profiles between low and high $I$-value. Because the reduced equations evolve only the phase differences, they permit denser parameter sweeps at substantially lower computational cost than simulations of the full delayed network. This efficiency makes it possible to resolve narrow stability regions, small-basin attractors, and bifurcation structures that can be difficult to identify using return maps alone.

Together, return-map analysis and phase reduction provide a unified framework for characterizing multistability, bifurcations, and transitions among rhythms in delayed inhibitory CPG networks. Their close agreement validates the phase-reduced description in the weak-coupling regime. The phase model should not, however, be viewed as a replacement for the full-system analysis: return maps remain essential when stronger coupling or delay-induced waveform changes produce amplitude effects that lie outside the assumptions of first-order phase reduction.

Future work should extend this framework to larger and less symmetric network architectures to determine how the number, stability, and organization of coexisting rhythms scale with network size and connectivity. A more detailed analysis of the local and global bifurcations responsible for the creation, stability changes, and disappearance of invariant curves would also clarify the mechanisms underlying transitions between phase locking, modulated rhythms, and persistent phase slipping.

ACKNOWLEDGMENT

X.Q.’s work partly funded by the National Natural Science Foundation of China (Grant Nos. 11872183) and China Scholarship Council (Grant Nos. 202506150077). A.L.S.’s research was partially funded by a National Science Foundation (NSF) (Award No. DMS-2407999).

DATA AVAILABILITY

Some of the code used for the numerical simulations is publicly available at [https://github.com/xinxin77-cloud/FHN3_phase_difference].

References