Neuromorphic design of a rhythmic network

Introduction

Rhythmic motor behaviours, including locomotion, respiration, and swallowing, rely on the repeated and coordinated activation of muscles. Biological systems typically govern these movements through central pattern generators (CPGs): neural circuits capable of producing rhythmic motor commands without requiring continuous rhythmic input [2, 3]. While sensory feedback and descending signals modulate the resulting behaviour, CPGs provide the fundamental organizing principle for motor control [4, 5]. Classical interpretations of biological CPGs distinguish between rhythms driven by intrinsic pacemakers and those emerging from network interactions [6, 7]. Pacemaker-driven architectures build upon intrinsically oscillatory neurons coordinated through synaptic coupling [8]. Conversely, network-driven circuits do not require individual neurons to oscillate in isolation; the rhythm arises entirely from structural interactions like reciprocal inhibition and post-inhibitory rebound [7]. Actual biological networks likely blend these mechanisms, with individual neuronal contributions shifting based on operating points, synaptic inputs, and neuromodulatory states [8]. Such architectures have deeply influenced robotic locomotion by providing a low-dimensional parametrisation for coordinated, high-dimensional movement [9, 10].

Traditionally, engineers have favoured the design of rhythmic actuator trajectories from the theory of oscillator networks [11, 12]. Translating the biological CPG concept into robotics predominantly yields networks of coupled limit-cycle oscillators, whether formulated through abstract models (e.g. Hopf, Van der Pol) or paired half-centers abstracted as effective oscillatory units (e.g., Matsuoka circuits) [13, 14]. In these architectures, the continuous nodal dynamics and coupling topology directly dictate synchronization patterns, phase relationships, and viable gait transitions. However, endogenous rhythm generation, while crucial, solves only one piece of the practical locomotion puzzle [15, 16, 17]. Movement initiation, adaptation, and termination critically rely on sensory feedback, reflexes, body mechanics, and descending commands [18, 19].

Consequently, robotic and neuromechanical designs either investigate purely reflex-driven architectures [20, 21] or augment feedforward CPG signals with proprioceptive feedback loops [22, 23, 24, 25, 26], exteroceptive sensing [27], and local reactive reflex pathways [28, 29, 30, 31]. Yet, the ideal method for integrating internally generated rhythms with incoming sensory data remains an open design challenge, particularly in asynchronous or event-based control schemes.

The event-based approach of neuromorphic control provides an alternative angle of attack. The present paper revisits CPG design from the perspective of interconnecting excitable rather than oscillatory nodes, focusing on tuning and robustness while bypassing the physiological details of biological neurons. A loop of nodes that propagate events forms a feedback system that does not separate the mechanism of endogenous rhythm generation from the interaction with the environment. Event-based design is hierarchical, and it includes automation (i.e., rhythm generation) and regulation (i.e. interaction with the environment) at any scale [32, 33, 34, 35, 36]. Formulating rhythmic motor control through discrete events structurally unifies the question of endogenous rhythm generation and the question of exogenous interaction, allowing nodal interactions and external perturbations to share the exact same event-based representation.

Excitable neurons are natural elementary units for such event-propagating networks [37]. Unlike an oscillator, the isolated behaviour of an excitable neuron is a stable equilibrium, the resting state, until an input drives the behaviour above threshold, thereby triggering a transient event such as a neuronal spike or burst [38, 39]. While the same neuronal model can be excitable or oscillatory in distinct parameter regimes [40, 41], these two operating states enforce fundamentally different mechanisms of network coordination. Oscillatory nodes possess intrinsic frequencies prone to incompatibility and conflict under parameter heterogeneity. As a consequence, networks of oscillators demand precise parameter balancing, forcing the joint tuning of nodal and network parameters to achieve stable phase-locking. Excitable nodes, in contrast, lack an autonomous rhythm. As a consequence, they introduce no competing intrinsic frequencies that require global reconciliation. Global rhythmic behaviour instead emerges purely from the composition of local discrete-event propagation, as each node independently processes incoming events into outgoing responses [42]. This isolates the tuning process, replacing the precise balancing of high-dimensional coupled differential equations with strictly local, independent design rules.

Post-inhibitory rebound provides a highly effective mechanism for organizing this sequential event propagation [43, 44, 45]. Because release from inhibition can actively trigger a spike or burst, inhibitory synapses do more than merely suppress activity; they structure alternation and drive sequential propagation. Specifically, directed ring topologies of cyclic inhibitory connections can sustain robust rhythmic patterns even when every constituent neuron remains entirely quiescent in isolation. Network and population models demonstrate that sequences, travelling waves, and low-dimensional rotational dynamics—from motor gaits to respiratory rhythms [46]—readily materialise without intrinsic pacemakers [47, 48, 49]. While many of these existing models capture collective activity through continuous firing-rate dynamics, our approach preserves the discrete, spike-based nature of neuronal signalling to frame rhythmicity purely as event propagation between excitable nodes. We mirror the concept of wave propagation in excitable media [50], but channel it through discrete inhibitory interactions and post-inhibitory rebound within an explicitly defined, discrete network topology.

Traditionally, the analysis tools of rhythmic networks have relied on the theory of limit-cycle oscillators. In control theory, describing function analysis and harmonic-balance methods have been central to predict oscillations in feedback loops [51, 52, 53, 54, 55, 56]. Entrainment and synchronization questions have been instead studied via the phase-response curve analysis [57]. Phase reduction maps an autonomous limit-cycle oscillator onto a single phase variable, and phase response curves (PRCs) quantify how external perturbations advance or delay the underlying cycle [58, 59, 60, 61]. This approach has been extensively adapted for various oscillator models, including relaxation systems [62, 63], excitable media via isostables [64], and functional network responses [39, 65]. By isolating phase interactions, researchers can tractably predict synchronization, phase-locked states, and overall stability in coupled oscillatory networks [66, 67, 68, 69, 70, 71]. Yet, standard describing functions and phase response curves presuppose autonomous limit-cycle oscillations, making them inapplicable to excitable nodes that settle at a stable resting state. Sieling et al. [72] first demonstrated that phase response concepts could be extended beyond autonomous limit cycles to individual non-oscillatory neurons. Our recent work [1] shifted this perspective from isolated cells to fully excitable, non-oscillatory networks by introducing the event describing function (eDF) to quantify nodal propagation delay and the event phase response curve (ePRC) to quantify phase sensitivity. The present article develops this eDF and ePRC foundation into a systematic design and regulation methodology for heterogeneous CPG architectures. The rest of the article is organised as follows: § 2 uses a clock analogy to contrast coupled-oscillator and excitable networks, showing how event-propagating architectures decouple local timing adjustments and accommodate parameter heterogeneity. § 3 addresses endogenous rhythm generation by developing the event describing function (eDF) as a design tool. We demonstrate how nodal parameter variations shape the eDF, enabling the systematic design of heterogeneous networks with prescribed properties, such as the total network period and phase profile. § 4 studies exogenous rhythm regulation using the event phase response curve (ePRC). We show how nodal parameters shape the ePRC under excitatory and inhibitory perturbations, providing a tractable local methodology for designing sensory sensitivity and network entrainment. Finally, § 5 summarises the results and outlines directions for future work.

Motivating example

The core premise of this work lies in a shift from classical coupled-oscillator networks to excitable networks. This structural distinction is illustrated via a Hodgkin–Huxley neuron in Figure 1: panel (a) displays an autonomous oscillatory regime generating a continuous limit cycle, whereas panel (b) shows an excitable regime that remains quiescent at rest until an incoming input event triggers an output event. Mathematically, the two distinct operating regimes of the same model only differ by a constant input current. Yet, the input–output properties of the two systems differ considerably, and they impact both the analysis and the design of the neuron in a neural network. In particular, the excitable regime is a sweet spot between a system that can generate events autonomously and a system that interacts with its environment.

To suggest the practical differences between oscillator-based and network-based CPGs, we start from the elementary design question of interconnecting neurons to produce a clock-like cyclic rhythm, which will serve as a guiding analogy for the remainder of the paper. We first demonstrate that decoupling local from global behaviour grants the excitable network superior robustness against heterogeneity. Following this, we contrast both paradigms directly by illustrating how a localised synaptic modulation propagates globally across the oscillatory network, yet remains isolated to local dynamics in the excitable ring.

Network design of a clock: sequential ring activation

In the absence of interaction with the environment, a clock is a dynamical system with a stable limit cycle characterised by the sequential activation of nodes in a predetermined firing order. In an excitable network, a unidirectional ring topology enables this objective with minimal connectivity (Figure 2, left). When an active node fires, it excites its downstream neighbour, establishing an event-driven propagation motif along the ring. In this framework, a nodal spike serves as a discrete trigger for the state transition of a discrete automaton. The continuous duration of an active state is governed by the delay between successive spikes. We can regard this motif as an abstraction of many biological CPGs, where neurons burst one after the other, and where the duration of the burst physically controls the duration and delay of the successive nodal events. We will return to the extension of the neuromorphic framework to bursting networks in the conclusion.

The left panel of Figure 2 illustrates the directed ring topology. We adopt a four-node network ($N = 4$) throughout this paper to maintain clear and legible visualizations. This architecture can be interpreted as clustering the conventional 12-hour clock positions into four functional groups, an aspect detailed in Section 3.4.3. The design principles presented here extend directly to rings of any size $N$.

The right panel of Figure 2 shows the resulting voltage traces, highlighting the activity of node $n_0$. As illustrated, these excitable nodes rest at a stable subthreshold equilibrium until an incoming inhibitory event hyperpolarises them, triggering a post-inhibitory rebound spike upon release (this mechanism is explained in detail in § 3.1. In this work, we focus on this fundamental ring topology, as it captures the fundamental event-driven foundation onto which more complex layers can be added.

Attempting to realise a ring oscillator with oscillatory nodes reveals fundamental limitations. When coupled with purely excitatory synapses, autonomous oscillators tend to lock in phase, collapsing the desired sequence into global synchrony. Conversely, inhibitory coupling can support phase-lagged sequential firing, but this solution is rarely unique; inhibitory oscillator rings frequently exhibit complex multi-stability [12]. Consequently, the emergence of the correct sequence becomes heavily reliant on specific initial conditions, leaving the desired rhythm vulnerable to parameter variations and transient perturbations.

To eliminate competing attractors and enforce the sequential firing order as the unique stable limit cycle in an oscillatory network, additional structural constraints are required. Specifically, one must combine global winner-take-all (WTA) inhibition—to ensure solitary activation—with a unidirectional excitatory ring—to enforce the firing sequence. An analogue of this architecture is seen in rebound winner-take-all (RWTA) networks [32], which apply this dual-layer topology to rebound-excitable neurons. This composite architecture is straightforward to implement with excitable neurons, and can yield robustness benefits relative to the ring topology. Yet, a similar oscillatory network remains inherently fragile, requiring global and precise tuning of intrinsic frequencies and synaptic gains to maintain coherent phase-locking [73]. As we will discuss in § 2.2, this makes the design problem similarly vulnerable to parameter variations and perturbations.

Robustness to heterogeneity

Robustness to parameter variation is essential for a straightforward design that avoids precise fine-tuning. Beyond simplifying implementation, this robustness directly dictates the network’s tunability: the system’s behaviour must remain qualitatively stable under parameter variations to allow for reliable modulation and control. With this in mind, let us first examine how the two architectures compare in terms of robustness to parameter heterogeneity. We deliberately benchmark the minimal excitable network against the most robust oscillatory configuration capable of producing a stable clock rhythm.

To quantify this robustness, we apply relative Gaussian parameter variations $p_i \sim\mathcal{N}(\bar{p}_i,(h \cdot\bar{p}_i)^2)$ to key network parameters ($g_{\mathrm{Na}}$, $g_{\mathrm{K}}$, and $\tau_{\mathrm{syn}}$). We simulate 400 trials for each heterogeneity level $h$, classifying a run as successful only if it yields a sustained, sequential pattern in the predetermined firing order. Figure 3 shows the fraction of these successful trials for both architectures. It demonstrates that, although the underlying neuronal and synaptic dynamics are (nearly) identical for the two architectures, the oscillatory network quickly exhibits pronounced fragility for heterogeneity levels above $h > 10\%$. On the other hand, the excitable network maintains its functional rhythm across a much wider range of $h$ before its validity eventually deteriorates.

This highlights a core advantage of the event-driven approach: timing differences remain localised and do not interact across the network, since pattern generation relies on threshold-based event propagation rather than continuous phase-matching.

Local modulation

We now investigate the impact of local parameter modifications on the resulting global activity in both architectures. While the spatio-temporal dynamics in each case are governed by the total network period $T_{\mathrm{tot}}$ and the internal phase profile $\Phi$, the mechanisms translating local changes into global modulation differ fundamentally.

In an excitable network, activity propagates strictly through discrete nodal delays $\Delta_i$. The overall network period is simply the sum of these delays, $T_{\mathrm{tot}} = \sum_{i=0}^{N-1} \Delta_i$, from which the internal phase profile $\Phi= \{\varphi_i\}_{i=0}^{N-1}$ follows directly as the set of relative delays, defined for each node by $\varphi_i = \Delta_i/T_{\mathrm{tot}}$. In an oscillatory network, however, each node possesses an intrinsic natural frequency, making the emergent global rhythm a complex product of coupled interactions.

Figure 4 demonstrates this crucial difference through a simple experiment, contrasting spike raster plots (Figure 4(a), (Figure 4b)) and effective inter-neuron delays (Figure 4(c), (Figure 4d)): at $t = 2000\ \mathrm{ms}$, we modify the network’s local behaviour by increasing the synaptic time constant $\tau_{\mathrm{syn}}$ from node $n_0$ to node $n_1$ in order to lengthen the delay between them from 20 ms to 40 ms. Within the clock rhythm framework, this operation directly corresponds to enlarging one specific interval without intending to alter the remaining sequence.

In the excitable network (Figure 4(a), (Figure 4c)), this modulation predictably alters only the local delay between node $n_0$ and node $n_1$. The remaining nodal dynamics are unaffected, proving that excitable networks provide decoupled local control over timing. In contrast, in the oscillatory network (Figure 4(b), (Figure 4d)), altering this single connection breaks the fragile phase-locking balance, causing unintended shifts in the effective delays across all nodes. Remember that, as discussed above, large levels of modulation may quickly result in the destruction of the rhythm in the oscillatory network.

Note that in (certain) coupled oscillator networks, a desired phase profile tuning can be achieved by accounting for the global effects of the synaptic connections [74]. However, it is the locality of the tuning in the excitable network that we want to stress here. Such local tuning is desirable in more complex networks or more complex modulation objectives.

Local analysis

These findings highlight why excitable networks offer a compelling alternative to networks of coupled oscillators: robust global coordination emerges from strictly local, nodal input–output characteristics. Because rhythm generation relies on decoupled event propagation, the global network dynamics collapse into a sequence of isolated, local transitions. This locality eliminates the need to analyse complex, high-dimensional coupled systems; instead, global network behaviour can be rigorously analysed, predicted, and designed purely through the input–output response of individual nodes. This principle of local analysis forms the primary theme for the remainder of this paper, establishing a fully distributed framework to analyse, design and tune event-driven rhythms.

Endogenous rhythm generation

In this section, we address the design problem for generating a rhythm in isolation, that is, in a network that does not interact with its environment. We postpone the interaction properties to the next section.

Mechanism of rhythm generation

We consider the design of rhythmic patterns in a network of excitable neurons interconnected by synapses. Excitable neurons have the ability to fire, that is, emit a voltage spike event, when triggered by a synaptic current (which itself is elicited by a presynaptic spike event). This is the core mechanism of event propagation in the network. If the presynaptic neuron is excitatory, the synaptic current depolarises the postsynaptic neuron, eliciting a spike if the depolarization is above threshold. If the presynaptic neuron is inhibitory, the synaptic current hyperpolarises the postsynaptic neuron, possibly eliciting an event when the inhibition is released. The latter phenomenon is called rebound excitability [75]. In the present paper, the primary focus is on event propagation in a network of inhibitory rebound neurons, as rebound propagation provides a robust mechanism to propagate events with adjustable delays, see e.g. [32].

Beyond the two classical currents of a mixed feedback neuron (fast activation of an inward current and slow activation of an outward current), two additional specific bio-inspired internal currents are used to control specific properties of the neuron. First, the rebound property is decoupled from excitatory spike generation and is controlled by a low-threshold, fast activation, slow de-inactivation, inward current ($I_{\mathrm{Na,Reb}}$) which is inspired by the low-threshold calcium current ($I_{\mathrm{Ca,T}}$), whose role has been studied extensively in the literature [45, 75, 76, 77]. Secondly, the latency of the spike generation mechanism is controlled by a fast activation, slow inactivation, outward current ($I_{\mathrm{A\text{-}potassium}}$), inspired by the biological A-type potassium current [78].

The full modelling details of the nodal and synaptic dynamics are provided in Appendix § 6.

The event describing function

The event map and event describing function

In the elementary ring network of Figure 2, each node acts as an adjustable delay between the pre- and post-synaptic events. The event map that defines this delay is shown in Figure 5: it consists of the series interconnection of the input synapse, modelled as an event-to-analog (E/A) converter; the neuron, modelled as a physical continuous-time circuit with input current $I_{\mathrm{syn}}(t)$ and output voltage $V(t)$; and the event readout, modelled as an analog-to-event (A/E) converter from the continuous-time signal $V(t)$ to a discrete event.

Vector (PDF)

Assuming a periodic input of presynaptic events, the event describing function (eDF) characterises the nodal delay as a function of the input period $T_{\mathrm{in}}$ [1]. This terminology is reminiscent from the describing function [51, 54] of a continuous-time system, when analysing continuous-time limit cycle oscillations in feedback systems.

Consider a periodic input event sequence $\mathcal{E}^{\mathrm{in}} = \{t_k^{\mathrm{in}}\}_{k\in\mathbb{Z}}$ with period $T_{\mathrm{in}}$ (such that $t_{k+1}^{\mathrm{in}} = t_k^{\mathrm{in}} + T_{\mathrm{in}}$). The describing event function assumes that the node operates in a 1:1 phase-locked regime where each input event triggers exactly one output event (this is required for simple event propagation). In steady-state, the eDF maps $\mathcal{E}^{\mathrm{in}}$ to an output event sequence $\mathcal{E}^{\mathrm{out}} = \{t_k^{\mathrm{out}}\}_{k\in\mathbb{Z}}$ via a scalar delay $\Delta$, such that:

where $\Delta(T_{\mathrm{in}})$ characterises the steady-state absolute nodal delay as a function of the input period. To express this timing relationship as a relative phase shift within the oscillation cycle, we may also define a relative eDF as:

The eDF curve

The eDF curve $\Delta(T_{\mathrm{in}})$ is constructed by numerically integrating the nodal continuous-time system dynamics across a sweep of input periods $T_{\mathrm{in}}$, each time extracting the steady-state time shift $\Delta$. The resulting eDF curve reveals three distinct regimes, as shown in Figure 6(a):

For large input periods ($T_{\mathrm{in}} \ge T_r$), the nodal delay is the same as in response to a single isolated input event. This is due to the fading memory nature of the excitable neuron [79]: its memory of the past decays over time, meaning that only the last input event determines the output when the input events are sparse enough. For increasing input periods, the nodal delay converges to a constant baseline ($\Delta(T_{\mathrm{in}}) = \Delta_\infty$) and the relative phase shift decays asymptotically toward zero $(\lim_{T_{\mathrm{in}}\to\infty}\varphi(T_{\mathrm{in}})=0)$. For small input periods, the 1:1 entrainment fails. This can have various underlying causes: driving the system faster than $T_{\min}$ causes excitatory nodes to exhibit phase slipping or higher-order $N:M$ locking (firing $N$ times for every $M$ input cycles), while inhibitory nodes enter a state of sustained inhibition that completely suppresses spike generation. In the intermediate regime where $T_{\min}<T<T_{\mathrm{r}}$, the steady-state delay decreases to the asymptotic value $\Delta_{\infty}$ and may exhibit nonlinear effects.

The existence of $T_{\mathrm{r}}$ and the asymptotic constant regime is an important consequence of nodal excitability. In practice, this allows the eDF to be computed from a single input event provided the global network period is large enough (i.e., $T_{\mathrm{tot}}>T_{\mathrm{r}}$). By ensuring the network operates in this fully-reset regime, the local nodal delay decouples entirely from the global period, remaining fixed at $\Delta_{\infty}$. This yields a major computational advantage: instead of integrating the system dynamics over multiple cycles to reach steady state, the exact steady-state delay can be extracted from one isolated transient response.

At the network level, determining whether an autonomous oscillation can emerge simply becomes an algebraic question. For instance, for a ring network, a global period $T_{\mathrm{tot}}>0$ must be a solution of the network period equation:

From a design perspective, inhibition is the preferred mechanism of transmission because it provides a longer, highly controllable nodal delay with two distinct dimensions of tuning. As shown in Figure 6(b), the eDF behaves fundamentally differently depending on the synapse type. Excitatory delays are inherently brief, restricted by the fast time constant of the spike upstroke and only the latency current ($I_{\mathrm{A}}$) as control parameter. In contrast, inhibitory delays are governed by both the synaptic decay time ($\tau_{\mathrm{syn}}$) and the neuron’s internal rebound latency ($I_{\mathrm{A}}$). This dual-control mechanism yields a broad tuning range, enabling even small networks to easily satisfy the global period condition.

Tuning the event map

The local analysis of the eDF establishes a principled foundation for network design. To enable a design framework through local nodal tuning, it is enough to identify key parameters that control the nodal delay. We make a distinction between internal and external tuning: internal tuning adjusts the rebound latency of the excitable neuron ($I_{\mathrm{A}}$) whereas external tuning adjusts the synaptic time constant ($\tau_{\mathrm{syn}}$).

As introduced earlier, we exploit the well-known functionality of the A-type potassium current ($I_{\mathrm{A}}$) in order to tune the nodal delay internally. The slow inactivation of this current provides a direct control of the rebound latency (that is, the transient behaviour after release of inhibition). Tuning its maximal conductance ($g_{\mathrm{A}}$) affords direct, exponential control over the event delay without compromising overall excitability.

Figure 7 illustrates how $\tau_{\mathrm{syn}}$ and $g_{\mathrm{A}}$ provide complementary control parameters to shape the eDF: the external, synaptic time constant $\tau_{\mathrm{syn}}$ offers fine, linear tuning (panel (a)), while the internal, cellular A-type conductance $g_{\mathrm{A}}$ offers high-dynamic-range exponential tuning (panel (b)). This design procedure offers a simple illustration of how the neurophysiology of excitability can be leveraged and serve as a biological inspiration for neuromorphic control.

Local tuning of network properties

Tuning the global period

As a first illustration, we show how nodal parameters (the synaptic decay time constant $\tau_{\mathrm{syn}}$ and intrinsic $I_{\mathrm{A}}$-conductance $g_{\mathrm{A}}$) can be exploited to tune the global network period $T_{\mathrm{tot}}$. In prior work, we showed that the collective network period can be predicted by summing individual nodal delays, and how $g_{\mathrm{syn}}$ and $\tau_{\mathrm{syn}}$ shape the empirical eDF curve [1]. Extending this principle to global network control, we evaluate two regimes: homogeneous modulation across all $N$ nodes, and localised heterogeneous modulation, where a single target node is adjusted while the remaining nodes remain fixed.

As shown in Figure 8, both homogeneous and heterogeneous adjustments of the synaptic time constant $\tau_{\mathrm{syn}}$ and the intrinsic A-type conductance $g_{\mathrm{A}}$ serve as effective control knobs for frequency regulation. Specifically, Figure 8(a) demonstrates that scaling $\tau_{\mathrm{syn}}$ (both globally across all nodes and locally for a single target node) tunes the global network period $T_{\mathrm{tot}}$ in a smooth, nearly linear fashion. In contrast, Figure 8(b) shows that varying the A-type potassium conductance $g_{\mathrm{A}}$ governs $T_{\mathrm{tot}}$ across a much wider dynamic range, exhibiting an exponential scaling behaviour.

Furthermore, the transitions between nominal and modulated rhythms occur cleanly without unexpected transient dynamics, as shown in Figure 9. Specifically, Figure 9(a) and (Figure 9b) illustrate the smooth temporal evolution of inter-neuron delays under homogeneous and heterogeneous parameter updates, respectively (changing $\tau_{\mathrm{syn}}$ from 2.15 to 6 at $t = 700$). At the same time, the corresponding membrane potential traces in Figure 9(c) and (Figure 9d) confirm that spike generation remains phase-locked and stable throughout the transition, avoiding skipped cycles or bursting anomalies. This rapid, stable adaptation stems directly from the event-propagation architecture, establishing $\tau_{\mathrm{syn}}$ and $g_{\mathrm{A}}$ as robust control knobs for network period regulation.

Tuning the phase profile

A complete description of the network rhythm requires both its global period and its phase profile $\Phi$ (the set of relative nodal delays $\varphi_i$). Because these delays operate independently, modulating individual event delays $\Delta_i$ allows us to seamlessly reshape the duty cycle without destabilizing rhythmic event generation.

Because homogeneous nodes cannot generate asymmetric phase profiles, we consider heterogeneous tuning to achieve a desired phase profile $\Phi^*$ (e.g., representing a desired locomotor gait). As an example, we study the transition from a uniform splayed walking gait with phase offsets $\Phi_{\mathrm{walk}}^{\mathrm{offsets}} = \{0, 0.25, 0.5, 0.75\}$ to a galloping gait with offsets $\Phi_{\mathrm{gallop}}^{\mathrm{offsets}} = \{0, 0.4, 0.5, 0.9\}$ in the case of a quadrupedal CPG ($N = 4$ nodes corresponding to the four limbs $LF$, $LH$, $RF$, and $RH$). Assuming an underlying phase–task mapping handles node reassignment during the gait transition to establish the correct sequential ordering, the remaining transition is achieved by reshaping the phase profile from $\Phi_{\mathrm{walk}} = \{0.25, 0.25, 0.25, 0.25\}$ to $\Phi_{\mathrm{gallop}} = \{0.1, 0.4, 0.1, 0.4\}$.

Starting from the nominal walking gait ($a$) where $\Delta_i^a = 40$ ms, we consider two approaches to obtain the target galloping gait ($b$).

In the first approach, we preserve the global network period ($T_{\mathrm{tot}}^a = T_{\mathrm{tot}}^b = T_{\mathrm{tot}}$) and derive the target absolute nodal delays directly from the phase profile: $\Delta_i^b = \Phi_i^* T_{\mathrm{tot}}$. The result is shown in Figure 10(a), (Figure 10c), (Figure 10e), where $T_{\mathrm{tot}}$ remains fixed while the gait rhythm is modulated at $t = 2000$ (by dropping $g_A$ from 1.2 to 0.0 for nodes $n_0$ and $n_2$, increasing $g_A$ to 1.45 for nodes $n_1$ and $n_3$, and increasing $\tau_{\mathrm{syn}}$ from 2.15 to 2.4 for synapses $3 \to0$ and $1 \to2$) from $\Phi_{\mathrm{walk}}$ to $\Phi_{\mathrm{gallop}}$ (depicting the spike raster plot, absolute nodal delays $\Delta_i$, and relative nodal delays $\varphi_i$, respectively).

In the second approach, we keep the absolute delays of nodes $n_1$ and $n_3$ unaltered while modulating nodes $n_0$ and $n_2$ to obtain the target phase profile $\Phi^*$. Using the constraint $\Delta_i^b = \Delta_i^a$ for $i \in\{1,3\}$, the new total period becomes $T_{\mathrm{tot}}^b = \sum_{i=0}^{N-1} \Delta_i^b = \Delta_0^b + \Delta_2^b + \Delta_1^a + \Delta_3^a$. Substituting the phase relationship $\Delta_i^b = T_{\mathrm{tot}}^b \Phi_i^*$ for $i \in\{0,2\}$ gives $T_{\mathrm{tot}}^b = T_{\mathrm{tot}}^b(\Phi_0^* + \Phi_2^*) + \Delta_1^a + \Delta_3^a$. Finally, we solve for the updated nodal delays:

The second strategy is illustrated in Figure 10(b), (Figure 10d), (Figure 10f), where the delays of nodes $n_0$ and $n_2$ are shortened relative to the total period to achieve the desired profile $\Phi^*$ (implemented at $t = 2000$ by dropping $g_A$ from 1.2 to 0.0 for nodes $n_0$ and $n_2$, and decreasing $\tau_{\mathrm{syn}}$ from 2.5 to 0.5 for synapses $3 \to0$ and $1 \to2$).

These two example highlight how a few local parameters enable a simple and intuitive tuning of the network rhythm properties.

Modular design and phase clustering

A basic ring network naturally enforces a strict, non-overlapping sequence of events. To encode more sophisticated rhythmic behaviours, we can expand the spatial resolution of the rhythm simply by increasing the number of nodes $N$. We illustrate the role of this higher-dimensional architecture using the quadrupedal locomotor gait from the previous section.

As before, the simplest realisation is a 1:1 phase–task mapping using $N = 4$ nodes, where each node is assigned to one of the four limbs as before ($LF$, $LH$, $RF$, and $RH$). In this setup, a spike from a node triggers the onset of that limb’s active phase (e.g., the stance phase), and the delay to the subsequent node dictates the timing of the next limb’s activation. While simple, a strict 1:1 mapping implicitly enforces sequential, non-overlapping active states: the onset of one limb’s phase perfectly coincides with the termination of the previous one.

However, many locomotor gaits require muscle activations to overlap, or require the duration of a limb’s active phase to be tuned independently of the broader gait sequence. Because an excitable node outputs a single discrete event, a 1:1 mapping provides only one timing reference (e.g., task onset) per limb. To encode both the initiation and termination of a continuous action purely within the network’s timing dynamics, we can expand the network resolution by assigning a cluster of nodes to a single high-level task.

For example, one might make use of a 2:1 mapping, where two dedicated nodes govern a single task. The first node’s spike initiates the active state (e.g., muscle contraction), and the second node’s spike terminates it. By exposing the task duration as a tunable nodal delay, this clustering strategy decouples the duration of the activity from the delay to the subsequent task, granting independent control over the duty cycle of individual limbs.

If a gait additionally requires overlapping activations across multiple limbs, this can be achieved by defining the phase–task mapping such that a specific network phase explicitly represents a joint state (e.g., ‘both $LF$ and $RH$ are currently active’).

This modular way of clustering can be extended to arbitrary resolutions. For example, the continuous motion of the $LF$ limb could be decomposed into a sequence of distinct functional sub-phases: initiation (heel-strike), maximum load (mid-stance), and release (toe-off). By mapping a cluster of three nodes to this sequence, each sub-phase is exposed to local tuning. Consequently, a designer can selectively compress or stretch a specific sub-task (e.g., extending the stance duration to accommodate a heavier load) by tuning a single local delay, all while the event-driven architecture naturally preserves the integrity of the broader network rhythm.

Robustness and failure modes of the event map

The simplicity and robustness of the design framework provided by simple tools as the event Describing Function (eDF) rely on acknowledging specific assumptions. Most importantly, the assumption that the nodal input events are sparse enough ($T_{\mathrm{tot}} > T_{\mathrm{min}}$), such that the nodal delay is well-defined and easy to compute. The inhibitory nature of the network and the rebound property of the nodes are central to that assumption. The design of the nodal rebound dynamics (determined by the intracellular current $I_{\mathrm{Na,Reb}}$) is in that sense a critical ingredient for the robustness of the global network. Because the critical property is so localised, its robustness analysis is tractable, irrespective of the network size. The key to robustness and tunability is to design the intracellular currents in such a way that the nodal delay is sensitive to few of them (for tunability) and insensitive to most of them (for robustness). This design philosophy is a key inspiration from cellular neurophysiology.

Exogenous rhythm regulation

Practical rhythmic networks are rarely designed to operate in isolation. While they can generate a rhythm endogenously, they must also integrate sensory feedback from the external environment. Such exogenous interaction is critical for robust operation; for example, it allows a locomotor network to seamlessly adapt its phase to uneven terrain. This section introduces a second local analysis tool for this external analysis.

The event phase response

In the classical theory of networks of oscillators, the key analysis tool for external properties is the Phase Response Curve (PRC) [59, 60]. We introduce the event Phase Response Curve (ePRC) as the analog tool for a node that is excitable rather than oscillatory: starting from a nominal periodic regime, it measures the steady-state time shift induced by an additional perturbing periodic event sequence.

Definition

Starting from a nominal periodic regime (e.g. the outcome of an endogenous rhythm design such as in the previous section), we introduce an external input via an external input synapse. The underlying assumption is that a neural network interacts with its environment via synaptic currents, consistent with the dynamics of biological neural circuits. For the sake of illustration, we restrict the analysis to a nominal internal rhythm designed from an inhibitory network, as in the previous section.

Let $T_{\mathrm{nom}}$ denote the period of the nominal, unperturbed endogenous input event sequence driving the excitable node into a 1:1 phase-locked event oscillation. A periodic perturbing event sequence $\mathcal{E}^{p} = \{t_k^{p}\}_{k \in\mathbb{Z}}$ is introduced at a time offset $\delta_p \in[0,T_{\mathrm{nom}})$ within each nominal period through the external synapse, such that $t_k^{p} = t_k^{\mathrm{in}} + \delta_p$. We apply periodic rather than isolated perturbations to enable a stationary analysis, as in the previous section.

The output event time shift $\Pi(T_{\mathrm{in}},\delta_p)$ defines the event phase response:

where $t_{\mathrm{nom}}^{\mathrm{out}}$ and $t_p^{\mathrm{out}}$ represent the nominal and perturbation-shifted output event timings, respectively.

When $T_{\mathrm{tot}}$ is sufficiently large such that the system fully relaxes to its resting state between each input event ($T_{\mathrm{tot}} > T_r$), the ePRC can be safely and efficiently evaluated from a single event transmission. Furthermore, in this regime, perturbations arriving after the output event produce no phase shift. We therefore restrict the evaluation window to perturbations preceding the nominal input event, such that $\delta_p \in[-T_{\mathrm{nom}},\Delta_{\mathrm{nom}})$.

Evaluating $\Pi(T_{\mathrm{in}},\delta_p)$ across a sweep of time offsets $\delta_p$ constructs the complete event phase response curve for a given endogenous behaviour. In parallel to classical PRC analysis, this curve defines the system’s capacity for entrainment under frequency detuning ($T_p \ne T_{\mathrm{nom}}$). This relationship relies on three principles: first, to remain synchronized, the node must shift its timing by the exact period difference $\Delta T = T_p - T_{\mathrm{nom}}$ every cycle. Second, vertically shifting the ePRC by this $\Delta T$ allows us to find the stable phase-locking offsets at the new zero-crossings. Finally, the amplitude of the ePRC establishes the absolute entrainment bounds, as any detuning $\Delta T$ exceeding the curve’s peak will prevent phase-locking entirely.

Contrast to phase response curves analysis

To place the ePRC within the established literature, we briefly contrast it with related phase reduction theories.

Unlike the classical infinitesimal Phase Response Curve (iPRC), which evaluates direct state-trajectory pulses on autonomous limit-cycle oscillators, the ePRC characterises non-oscillatory, excitable systems driven by explicit synaptic dynamics.

To account for such physiological synapses, the functional Phase Response Curve (fPRC) utilises finite-amplitude, repetitive perturbations [65]. While the fPRC can successfully characterise non-oscillatory, rebound-excitable components, it operates in a closed-loop fashion, adaptively applying perturbations at a fixed delay after each output spike. In contrast, the ePRC evaluates responses to a strictly independent periodic sequence, yielding a purely open-loop, input–output mapping of forced phase-locking.

Finally, the infinitesimal Isostable Response Curve (iIRC) explicitly addresses the strongly relaxing, return-to-rest dynamics of excitable nodes that possess a stable fixed point rather than a periodic orbit [64]. While both the ePRC and the iIRC account for relaxation, the iIRC quantifies continuous-time trajectory shifts following an isolated state perturbation. The ePRC abstracts away the continuous return trajectory entirely, providing a purely discrete input–output timing map for driven event propagation.

Excitatory vs. inhibitory ePRC

The ePRC profile differs qualitatively between excitatory and inhibitory perturbations (with nominal coupling remaining inhibitory throughout). Most importantly, excitatory inputs consistently produce phase advance ($\Pi(T_{\mathrm{in}},\delta_p) < 0$), whereas inhibitory coupling consistently induces phase delay ($\Pi(T_{\mathrm{in}},\delta_p) > 0$): by hyperpolarising the membrane and prolonging the time needed to reach threshold, it delays the subsequent output event. This is a striking property of the event phase response curves in Figure 11. Note that this should not be interpreted as a mathematical property, but rather as a phenomenological consequence of the modelling assumptions: for the neuron and synapse model in this paper (see Appendix § 6), and assuming the perturbation does not perturb the steady-state 1:1 transmission behaviour of the neuron.

The stereotypical profile of the ePRC is shown for excitatory and inhibitory perturbations in Figure 11(a) and (Figure 11b), along with a sample simulation showing the effect of the perturbation on voltage trace in Figure 11(c) and (Figure 11d). The clean complementarity between phase advance and phase delay enables bidirectional tuning of network phase advances and delays. This is a considerable simplification and novelty with respect to the classical framework of oscillator network analysis.

We stress that specific neuronal models might exhibit regions of phase advance (delay) in their inhibitory (excitatory) ePRC, for specific choices of the parameters. The classical literature often associates Type 1 PRCs with strictly positive phase advances [73, 80, 81]. The strict absence of zero-crossings in both our excitatory and inhibitory ePRCs suggests that our underlying model exhibits Type 1 excitability [73]. In any case, our main message is to stress the value of this property from a design perspective, rather than analysing whether this property may fail in specific parameter regimes. With the models of this article, we observe this property as the rule rather than the exception, without fine-tuning of the parameters.

Phase resetting and entrainment dynamics

The ePRC governs the phase-resetting mechanism and dictates the stability of entrainment under external periodic forcing. The related phase transition map reveals critical phase transition points (e.g., borders of 1:1 phase-locking, regions of event skip, or multi-spike regimes).

We can study simple entrainment ($T_{\mathrm{nom}} = T_{\mathrm{p}}$) directly through the ePRC, following the standard approach of phase response theory: zero-crossings identify the fixed points of the phase transition map, while their slopes determine stability. From the ePRC obtained above, we can predict the state of entrainment of the node, given the configuration of the nominal and perturbing synapse and the respective intrinsic periods $T_{\mathrm{nom}}$ and $T_{\mathrm{p}}$. However, since neither ePRC possesses a zero-crossing, stable-phase-locking is not predicted. On the other hand, phase-locking is achieved at fixed points $\delta_{\mathrm{p}}^{*}$ satisfying ePRC($\delta_{\mathrm{p}}^{*}$) $= T_{\mathrm{p}} - T_{\mathrm{nom}}$ represents the period mismatch between the nominal oscillation and the perturbing signal.

Figure 12(a) shows the entraining effect of a perturbing input to the examined node ($n_3$) during network operation ($N = 4$) and a period mismatch $T_{\mathrm{p}} - T_{\mathrm{nom}} = -4\,\mathrm{ms}$ for an intrinsic network period $T_{\mathrm{tot}} = 60.9$. The stable zero-crossing of the excitatory ePRC, shifted by $T_{\mathrm{p}} - T_{\mathrm{nom}}$, is located at $\delta_{\mathrm{p}} \approx9$. As predicted, we observe stable phase-locking of the node to the perturbations, with a phase offset of $\delta_{\mathrm{p}}^{*} = 9.13$.

This entrainment experiment can be repeated over a range of $T_{\mathrm{p}} - T_{\mathrm{nom}}$ and $g_{\mathrm{syn}}$ to find the boundaries of phase-locking. In Figure 12(b), the entrainment region of the network is shown for excitatory perturbations. As expected, the entrainment region matches the maximum values of the ePRC ($\approx5.5$ for $g_{\mathrm{syn}} = 0.05$).

Local tuning of the phase sensitivity

Analogous to the single-node event map analysis from the previous section, we examine how intrinsic parameter variations reshape nodal phase sensitivity. Because the maximum and minimum values of the ePRC directly determine the bounds of the entrainment region (and thus dictate the node’s input sensitivity), modulating the amplitude of this curve provides direct control over the system’s susceptibility to external forcing.

Single-parameter tuning

We first evaluate the impact of two key parameters on the amplitude and shape of the ePRC. Just as in the previous section, we focus on $\tau_{\mathrm{syn}}$ and $g_{\mathrm{A}}$ for the sake of illustration.

Figure 13 shows the excitatory ePRC across sweeps of $g_{\mathrm{A}}$ (a) and $\tau_{\mathrm{syn}}$ (b). Consistent with the nodal eDF analysis, both parameters directly govern the nominal onset delay, indicated by the rightmost boundary where the ePRC returns to zero. However, modulation of the external (edge) parameter $\tau_{\mathrm{syn}}$ of the nominal synapse quasi-linearly delays the nominal output (as established in § 3) while simultaneously deepening the ePRC and shifting its peak to later offsets. In contrast, increasing the internal cellular parameter $g_{\mathrm{A}}$ delays, deepens, and widens the ePRC in an exponential manner. Interestingly, its peak sensitivity remains centred at essentially the same phase offset, since the node becomes receptive as soon as the nominal synaptic current subsides. This qualitative distinction can be explained by the fact that the recovery of the neuron from inhibition is very sensitive during the phase of active $I_{\mathrm{A}}$. As a consequence, small perturbations have a strong effect during this specific current activation. This is in direct contrast to the effect of $\tau_{\mathrm{syn}}$, which effectively lengthens the time of full inhibition, during which the neuron is not sensitive at all.

Covaried tuning

Ideally, we wish to tune the sensitivity of a node independently of its endogenous behaviour by co-varying $g_{\mathrm{A}}$ and $\tau_{\mathrm{syn}}$. Because increasing $g_{\mathrm{A}}$ enhances sensitivity at the cost of lengthening the nominal onset, concurrently decreasing $\tau_{\mathrm{syn}}$ compensates for this shift, yielding an adjustable sensitivity level while maintaining identical nominal timing.

In order to counteract the difference between the linear effect of $\tau_{\mathrm{syn}}$ and the non-linear effect of $g_{\mathrm{A}}$, we employ a Hill-type sigmoid modulation for $g_{\mathrm{A}}$:

where $g_{\mathrm{A}}^{-}$ and $g_{\mathrm{A}}^{+}$ are the lower and upper parameter bounds, $g_{\mathrm{A}}^{\mathrm{mid}}$ represents the midpoint value at $x = 0.5$, $b$ locks this midpoint constraint, and the exponent $a$ controls the steepness to balance the initial delay and late recovery. Here, the shaping parameters $a$ and $g_{\mathrm{A}}^{\mathrm{mid}}$ depend on the nodal dynamics and are fitted empirically. Concurrently, $\tau_{\mathrm{syn}}$ decreases linearly along the sweep index $x$:

where $\tau_{\mathrm{syn}}^{-}$ and $\tau_{\mathrm{syn}}^{+}$ are the lower and upper parameter values.

This tuning, enabled by the A-type potassium current and synaptic decay time constant ($I_{\mathrm{A}}$ and $\tau_{\mathrm{syn}}$), is showcased in Figure 14. As illustrated, sweeping the covariation index $x$ from 0 to 1 progressively deepens the phase advance, increasing the node’s sensitivity to perturbations. As desired, the rightmost boundary of the ePRC remains firmly anchored at a constant phase offset throughout the sweep. This confirms that the covariation accurately preserves the nominal onset delay, successfully decoupling the node’s exogenous sensitivity from its endogenous rhythmic timing.

Local tuning of the network interaction properties

We now extend our framework from single-node sensitivity to network-level exogenous regulation, using parameter heterogeneity as an explicit design tool.

For the sake of illustration, we consider again the ring topology, and model interaction with the environment via additive synaptic inputs spread over all the neurons of the network. Rebound neurons are mostly sensitive to external disturbances at the end of an inhibitory presynaptic event, meaning that only local perturbations can advance or delay the network rhythm. In other words, the network sensitivity is not different from the neuron sensitivity.

However, coupling across the ring can introduce an additional timing cascade: when a perturbation shifts the output timing $t_i^{\mathrm{out}}$ of node $n_i$, the input arrival time $t_{i+1}^{\mathrm{in}}$ at the subsequent node $n_{i+1}$ is shifted by the identical amount. This consequently updates the effective perturbation phase offset $\delta_{\mathrm{p}}^{i+1}$ for all downstream nodes in the sequence.

Framed through the analogy of quadrupedal locomotion, where sensory feedback must be selectively interface with specific limbs or gait phases, we examine three strategies for shaping this network sensitivity: (1) global homogeneous tuning, (2) single-node heterogeneous tuning, and (3) clustered functional tuning.

Homogeneous tuning

We first consider global homogeneous tuning, where all $N$ nodes receive identical sensory input (e.g., representing a shared feedback signal across all four limbs) and undergo identical parameter modulation.

We first illustrate how the entrainment of the network can be modulated without changing the nominal rhythm.

Figure 15(a) illustrates this mechanism: applying parameter covariation (increasing $g_{\mathrm{A}}$ while concurrently decreasing $\tau_{\mathrm{syn}}$) uniformly across all nodes enables the network to entrain to periodic inputs that fail to entrain the nominal, insensitive circuit. The uniform modulation scales the amplitude of every neuronal ePRC equally, increasing the sensitivity of the network across the entire ring.

Because all nodes share identical sensitivity profiles, each position in the cycle is equally receptive to external driving. This spatial symmetry creates $N$ multi-stable, phase-shifted entrainment equilibria, allowing the network to lock to incoming periodic signals at any of the $N$ phases.

Single-node heterogeneous tuning

We next investigate how the modulation of a single node in an otherwise uniform network can modulate the network sensitivity.

Following the analogy of quadrupedal locomotion, we consider phase-dependent sensory entrainment driven by periodic ground-contact events (such as a heel strike). Rather than distributing sensitivity across all limbs, natural gaits often require sensory responsiveness to be restricted to a specific phase of the locomotor cycle. Starting from a globally insensitive network, rendering node $n_0$ (representing the touch-down phase of a reference limb, such as the left forelimb) sensitive allows the entire network to phase-lock to the external sensory signal exclusively through that phase. This single sensitive node acts as a primary phase-locking anchor, illustrating how parametric heterogeneity can dictate global entrainment. This provides a computational analogue to phase-dependent reflex gating and sensory resetting in biological locomotion, where proprioceptive or tactile feedback entrains or resets the gait rhythm only when received within a narrow, functionally appropriate time window.

While the nominal network completely ignores the periodic stimulus, Figure 15(b) shows how modulating the sensitivity of a single node ($n_0$) enables entrainment of the entire network. Upon making node $n_0$ sensitive, a single stable fixed point emerges in the network phase response exclusively at that phase. Unlike the $N$ stable states produced by homogeneous tuning, this single-node modification selectively locks the network to a unique phase state, leaving all other phases insensitive to sensory perturbation.

Clustered and functional tuning

As a final illustration, we consider the modulation of a cluster of nodes to enable modular, multi-phase entrainment that maps directly to the clustered phase–task mapping established in § 3. Grouping nodes into functional modules allows specific phases of an execution sequence to remain highly receptive to external feedback while keeping other phases autonomous and insensitive.

Returning to our quadrupedal locomotion analogy, the continuous muscular activity of a single limb can be partitioned into distinct sub-phases, each representing a different muscle activation profile during the stance phase. By assigning a dedicated sequence of nodes to these sub-phases, their sensitivities to external inputs can be modulated entirely independently. For instance, the initial impact sub-phase of the limb could be rendered insensitive to avoid disruption from transient mechanical noise upon touch-down, while a subsequent load-bearing sub-phase is tuned to be highly sensitive. This functionally gated entrainment window allows incoming environmental cues to reliably phase-lock the rhythm only when the limb is dynamically ready to respond.

We can also revisit the two-node functional clustering strategy introduced in Section 3.4.3. While the basic framework in this work treats tasks as point events (or as remaining active until superseded by the next event), certain motor commands require explicit, fixed execution windows. By assigning dedicated nodes within a clustered phase, we can decouple phase transitions from phase durations. Specifically, using two nodes per limb phase allows a (sensitive) ‘transition node’ to govern onset delay and entrainment bounds in response to external inputs (e.g., adapting touch-down timing to uneven terrain), while a downstream, (insensitive) ‘duration node’ dictates the active execution time of the task (e.g., maintaining a fixed muscle contraction duration). This maintains the simplicity of discrete event propagation rather than relying on continuous modulation of event durations.

Much like the endogenous design problem, the clustered sensitivity strategies highlight the modular nature of the architecture. They allow complex, multi-stage pattern generators to selectively interface with environmental signals only when sensory feedback is functionally relevant, without altering the underlying network connectivity.

Conclusion

In this work, we presented a distributed, event-based framework for rhythmic pattern generation in Central Pattern Generators (CPGs) composed of excitable, non-oscillatory neurons. Our approach rests on a conceptual shift, adopting the perspective that biological rhythms can emerge from sequences of discrete events in excitable neuronal networks rather than from coupling between neuronal oscillators. Our main goal was to illustrate that the the analysis and design of such rhythmic sequences is modular, robust, and scalable.

We addressed separately the internal question of generating a particular rhythmic sequence in a network and the external question of tuning how the internal rhythm can interact with the external environment. The internal question was analysed in § 3 by means of the event Describing Function (eDF). We showed how to modulate the global period and relative phases of a cyclic sequence through local control of the nodal transmission delays ($\Delta_i$). The external question was analysed in § 4 by means the event Phase Response Curve (ePRC). We showed how to modulate the entrainment properties of a cyclic network by local modulation of the nodal sensitivity. The results illustrate that event-driven excitable architectures naturally accommodate heterogeneity, enabling flexible phase–task mappings, functional node clustering, and robust rhythm generation.

All the illustrations in this paper were provided on the simplest example of rhythm, namely the generation of a periodic sequence in a ring network. This restriction was adopted for the sake of illustration, but the proposed methodology is general and should be investigated in more general architectures. In particular, two specific problems deserve further attention for the application of rhythm generation to rhythmic robotics. First, the analysis tools (eDF and ePRC) should be generalised to bursting (rather than spiking) neurons. Second, the generation of periodic sequences should be considered in network topologies that combine excitation and inhibition, such as the rebound winner-take-all networks studied in [32]. The results in this paper provide encouraging first steps for the general question of generating hierarchical event sequences in neuromorphic control systems.

Funding

The research leading to these results has received funding from the European Research Council under the Advanced ERC Grant Agreement SpikyControl n.101054323.

Declaration of Generative AI in the Writing Process

During the preparation of this work, the first author used Gemini (Google) to assist with language polishing and readability enhancements.

Data availability

The simulation code supporting this work are openly available on GitHub and permanently archived on Zenodo [82] (https://doi.org/10.5281/zenodo.22939323).

Mathematical models

Discrete-event node

The discrete-event nodes framing our rhythmic network are implemented as a three-stage functional pipeline: an event-to-analog (E/A) synapse, a continuous-time conductance-based neuron, and an analog-to-event (A/E) detector (as illustrated in Figure 5). When discrete event timings arrive at the node, the E/A synapse converts them into a continuous synaptic current. This current drives the excitable neuron’s state equations, generating a continuous membrane potential trajectory. Finally, the A/E event detector records a new discrete output event whenever this membrane potential crosses a predefined threshold.

Conductance-based mixed-feedback model

The core excitable node used throughout the network simulations (Figure 2–Figure 15) is a reduced two-dimensional conductance-based model governed by the mixed-feedback motif. The system state is defined by a fast membrane potential variable $V_f$ and a slower recovery variable $V_s$:

where $C_m$ is the membrane capacitance, $\tau_s$ is the slow time constant, $I_{\mathrm{app}}$ is the externally applied current, and $I_{\mathrm{syn}}$ is the sum of the incoming synaptic currents. The intrinsic ionic currents are formulated as:

Notice that we structurally decouple the transient sodium current into two distinct components: a high-threshold excitatory current ($I_{\mathrm{Na}}$) and a low-threshold rebound current ($I_{\mathrm{Na,Reb}}$). This allows for independent tuning of direct excitability and post-inhibitory rebound behaviour.

For simplicity of tuning and neuromorphic design, the voltage-dependent activation ($m_x$, $n_x$) and inactivation ($h_x$) gating variables are defined using a generalised sigmoid function $S(x)$. While any standard sigmoidal curve would suffice, we explicitly define $S(x)$ as a saturating rectified linear unit (SatReLU) to provide strictly linear control boundaries:

The gating variables are then constructed using independent bias ($V_{\mathrm{bias}}$) and slope ($k$) parameters. Some gates ($m_{\mathrm{Na}}$, $m_{\mathrm{A}}$ and $m_{\mathrm{Na,Reb}}$) depend on the fast voltage $V_f$, while others ($h_{\mathrm{Na}}$, $h_{\mathrm{A}}$, $h_{\mathrm{Na,Reb}}$ and $n_{\mathrm{K}}$) depend on the slow recovery voltage $V_s$:

Unless otherwise specified in specific experiments (e.g., during parameter sweeps or structural modulations), the nominal baseline parameters for the node are set as follows. The membrane parameters are $C_m=1.0\,\mu\mathrm{F}/\mathrm{cm}^2$ and $\tau_s=5.0\,\mathrm{ms}$. The maximal conductances (in $\mathrm{mS}/\mathrm{cm}^2$) are $g_{\mathrm{Na}}=50.0$, $g_{\mathrm{K}}=1.0$, and $g_{\mathrm{L}}=0.1$, with reversal potentials (in mV) of $E_{\mathrm{Na}}=100.0$, $E_{\mathrm{K}}=-20.0$, and $E_{\mathrm{L}}=0.0$. The A-type potassium current operates with $E_{\mathrm{A}}=-20.0$, while $g_{\mathrm{A}}$ acts as a primary tuning variable (typically $g_{\mathrm{A}}=0.0$ for the baseline node, and $>1.0$ when specifically tuning transmission delays).

The gating kinetics (bias in mV, slope in $\mathrm{mV}^{-1}$) are:

• Na: activation bias 10.0, slope 0.02; inactivation bias 30.0, slope 0.05.

• Na,Reb: activation bias $-10.0$, slope 0.1; inactivation bias $-5.0$, slope 0.1.

• K: activation bias 10.0, slope 0.05.

• A: activation bias $-15.0$, slope 0.1; inactivation bias $-10.0$, slope 0.15.

To transition the node from the excitable to the oscillatory regime, the rebound conductance is disabled ($g_{\mathrm{Na,Reb}} = 0.0$), the A-type conductance is elevated ($g_{\mathrm{A}} = 5.076$) to retain consistent intrinsic network period for $N = 4$, and a constant tonic drive of $I_{\mathrm{app}} = 4.0\ \mu\mathrm{A}/\mathrm{cm}^{2}$ is applied to bias the neuron into the oscillatory regime.

Hodgkin-Huxley neuron model

For the direct comparison of classical continuous dynamics versus event-driven dynamics in Figure 1, we employed a standard Hodgkin-Huxley model augmented with an A-type potassium current. The state equation is:

The gating variables evolve according to classical $\alpha/\beta$ exponential kinetics, which we shifted by $+65\ \mathrm{mV}$ to align the nominal resting potential near $0\ \mathrm{mV}$. The maximal conductances are $G_{\mathrm{Na}} = 120.0$, $G_{\mathrm{K}} = 36.0$, $G_{\mathrm{A}} = 10.0$, and $G_{\mathrm{L}} = 0.3\ \mathrm{mS}/\mathrm{cm}^{2}$, with reversal potentials $E_{\mathrm{Na}} = 115.0$, $E_{\mathrm{K}} = -12.0$, and $E_{\mathrm{L}} = 10.6\ \mathrm{mV}$. The autonomous oscillatory regime is driven by a constant tonic current of $I_{\mathrm{app}} = 15.0\ \mu\mathrm{A}/\mathrm{cm}^{2}$, whereas the excitable regime receives $I_{\mathrm{app}} = 0.0$ alongside discrete synaptic event triggers.

Event-based synapse

Incoming discrete events trigger a dynamic synaptic conductance. To ensure smooth integration while retaining switch-like activation, the synaptic drive is governed by a logistic sigmoid based on the presynaptic membrane potential $V_{\mathrm{pre}}$:

This drive activates the synaptic gating variable $s$ via a scaled rise and decay mechanism:

where $\alpha= \max\left(0,\left(\frac{1}{\tau_{\mathrm{rise}}}-\frac{1}{\tau_{\mathrm{syn}}}\right)d(V_{\mathrm{pre}})\right)$. The resulting applied current on the postsynaptic neuron is $I_{\mathrm{syn}} = g_{\mathrm{syn}}s(V_{\mathrm{f}} - E_{\mathrm{rev}})$.

By default, network interactions rely on inhibitory coupling with $g_{\mathrm{syn}} = 3.0\ \mathrm{mS}/\mathrm{cm}^{2}$, $E_{\mathrm{rev}} = -20.0\ \mathrm{mV}$, $V_{\mathrm{th}} = 30.0\ \mathrm{mV}$, a steepness of $k_{\mathrm{syn}} = 0.1$, a fast rise time $\tau_{\mathrm{rise}} = 0.1\ \mathrm{ms}$, and a baseline decay time $\tau_{\mathrm{syn}} = 2.15\ \mathrm{ms}$. Sensory perturbations (used in ePRC characterization) utilise a maximal synaptic conductance of $g_{\mathrm{syn}} = 0.05\ \mathrm{mS}/\mathrm{cm}^{2}$ with $E_{\mathrm{rev}} = 50.0\ \mathrm{mV}$ for excitatory inputs, and $g_{\mathrm{syn}} = 0.15\ \mathrm{mS}/\mathrm{cm}^{2}$ with $E_{\mathrm{rev}} = -20.0\ \mathrm{mV}$ for inhibitory inputs.

Event detection

The analog-to-event (A/E) conversion acts as the final stage of the nodal pipeline. A discrete output event is registered at time $t$ when the fast membrane potential $V_{\mathrm{f}}(t)$ crosses a detection threshold $V_{\mathrm{detect}} = 30.0\ \mathrm{mV}$ (or $50.0\ \mathrm{mV}$ depending on the amplitude of the chosen parameter regime) with a strictly positive derivative. These discrete timestamps formulate the outgoing event sequence $\mathcal{E}^{\mathrm{out}}$ that propagates to downstream nodes.

Simulation details

All neuronal models and network architectures were simulated computationally using custom Python scripts. To accommodate the specific numerical requirements of different analyses, integration was performed using a Forward Euler scheme. For extended network simulations (figures 2, 3, 4, 8, 9, 10, 12, and 15), a time step of $dt = 0.01\ \mathrm{ms}$ (or up to $0.025\ \mathrm{ms}$ where stability permitted) provided sufficient resolution. However, evaluating the discrete-event mapping curves (the eDF and ePRC in figures 6, 7, 11, 13, and 14) required capturing minute, sub-millisecond phase shifts. For these characterizations, the integration step was reduced to $dt = 0.001$ or $0.0001\ \mathrm{ms}$. Furthermore, to strictly eliminate discretization artifacts when registering event timings, the exact moment of threshold crossing ($V_{\mathrm{detect}}$) was calculated between integration steps using linear interpolation. Simulation routines were accelerated via Just-In-Time (JIT) compilation using the Numba library.