The neural circuit for stick insect locomotion produces rhythmic activity that drives locomotion in the intact animal as well as a fictive locomotor rhythm, which emerges when the circuit is isolated and hence feedback signals from stepping, ground contact, and so on are not present. Our first goal was to tune the parameters of the mesothoracic model to elicit an idealized time course of sequential unit activation that has been reported in the literature during the a single-leg locomotor rhythm when interactions between joint blocks associated with intact or semi-intact leg stepping are present (Büschges et al., 1994; Fischer et al., 2001; Akay et al., 2004; Daun et al., 2009; Bidaye et al., 2018); to maintain focus, we consider one specific rhythm, although a range of rhythms will certainly occur in practice (Graham and Wendler, 1981; Graham, 1985; Rosenbaum et al., 2010; Dallmann et al., 2019). The rhythm of interest features the following properties, which we henceforth call (P):

1.

The pattern is periodic in timeFootnote 1.

2.

Within each joint core, the two interneuron units take turns being active.

3.

The activation onsets of the Ext and Lev units are sufficiently synchronized.

4.

The Pro unit’s activation occurs approximately halfway through the active phase of Lev.

5.

The Lev unit’s activation ends approximately halfway through the active phase of Pro.

6.

The ends of the active phases of the Ext and Pro units (equivalently the activation onsets of Flx and Ret) are sufficiently synchronized.

We note that attaining these properties from the coupling architecture shown in Fig. 1 is non-trivial. For example, when Lev activates, it excites both Pro and Ext; these two targets must nonetheless activate at different times, while still deactivating together, to satisfy (P).

In the following, we first show results from the circuit depicted in Fig. 1 (right panel), where motoneurons are excluded, and the inter-joint excitatory couplings originate from an IN pool, as an extreme simplification of the feedback signaling involved in rhythmic leg stepping (cf. Section 2). These results include simulations obtained using XPPAUT (Ermentrout, 2002) and MATLAB and analysis based on the concepts and frameworks presented in Section 3. We explain the mechanisms underlying the rhythms generated by this circuit and discuss the robustness of its parameters. Subsequently, we extend the circuit by incorporating inter-joint inhibitory coupling in Section 4.4 and motoneurons in Section 4.5. We then evaluate how these extensions may influence the generation of a robust and ideal rhythm (P), identifying potential advantages and complications that they introduce.

Fig. 5

A baseline locomotor rhythm with timing properties (P). Top panel: Voltage time courses of INs, matching properties observed in stick insect stepping (Büschges et al., 1994; Fischer et al., 2001; Akay et al., 2004; Daun et al., 2009; Bidaye et al., 2018). The top set of traces show a set of INs that interact synergistically to achieve leg lift and advance (swing), while the bottom traces show their respective antagonists (for stance) in matching colors (Fig. 1). In the grey phase (2), the stance INs are active; in the darker green phase (4), the swing INs are active; and the pale green phases (1), (3) represent intermediate transitions. The model parameter values used are specified in the columns for “Escape" or “only IN, Escape" of Tables 1-5 in Appendix A. Bottom panel: The transition from active (bold color) to silent (pale color) phase or silent to active phase for each IN unit over one step cycle

4.1 A limb locomotor neural circuit model produces an escape-based stepping pattern

We consider a neural rhythm generated by a mesothoracic limb locomotor neural circuit composed of three coupled joint blocks, each including a joint core of two INs interacting by reciprocal inhibition. This rhythm features multiple phases of prolonged activation of different INs, and since the inactive IN in each block is subject to inhibition from its active partner, either escape or release could be associated with each phase transition in this rhythm. In practice, we found that producing a stable rhythm with the desired timing and phase relations required a tuning that featured phase transitions by escape. The activation time courses of the six IN units in the circuit for a baseline parameter set (see Appendix A), with transitions by escape, are shown in Fig. 5. Note that in an intact animal, activation of Lev drives levator muscles that lift the leg, subsequent activation of Pro drives protractor muscles that move the leg forward, and finally activation of Dep drives depressor muscles that return the leg to the ground; thus, a locomotor swing phase can roughly be defined as the time from the start of Lev activation until the end of Pro activation, which is approximately 40\(\%\) of the total cycle duration (Fischer et al., 2001; Grabowska et al., 2012). We next turn to phase plane analysis to explain how the features of this rhythm emerge and to lay the groundwork for understanding robustness of this rhythm to parameter variations.

4.2 Nullcline analysis identifies three classes of transitions within escape-based rhythms

Although the mesothoracic model comprises a system of 12 coupled nonlinear ODEs, a helpful alternative to considering a 12-dimensional phase space is to visualize its trajectories in a set of six phase planes, one per CPG unit, each including the projection of the trajectory along with multiple sets of v-nullclines corresponding to different levels of inputs that affect the visualized unit at different times during a locomotor cycle (Section 3.1; Fig. 6, with zoomed-in views shown in the outer columns). Specifically, each unit receives inhibitory input from its antagonistic partner within its joint core; the green pair of nullclines arise when inhibition is on and the blue (not visible in the zoomed-in plots) when it is off. Moreover, each unit receives excitatory input from one or more sources; in each same-color pair of nullclines, the solid nullcline arises when excitation is on and the dashed when it is off (e.g., green pair in zoomed view in Fig. 6). We will denote the attracting nullcline branches and associated structures for unit i, with \(i \in \\), by \(v^i_L(h), v^i_R(h), (v_^i,h_^i)\), and so on; in doing so, we drop explicit reference to the input levels on which these values depend, although we will occasionally include these as needed.

Fig. 6

The (v ,  h) phase planes for six coupled units, along with their corresponding zoomed versions. Two pairs of cubic-shaped v-nullclines are depicted: the green pair of nullclines arises when antagonistic inhibition is on and the blue pair when it is off. In each same-color pair, the solid nullcline represents the scenario when excitatory inputs are on, and the dashed nullcline indicates when they are off. The dotted gray curve denotes the h-nullcline. Each phase plane includes the projection of a limit cycle, colored according to the scheme in Fig. 5. A zoomed view near the left knees of the inhibited nullclines is included for each phase plane to provide a more detailed visualization. Blue labels indicate transition types

Each zoomed nullcline visualization in Fig. 6 focuses on the left knee, \((v_^i,h_^i)\), where the projected trajectory lies when the corresponding unit begins its transition from inactive to active, and hence reveals details of how each unit activates. For each projection, we refer to unit i as silent or inactive when its trajectory lies on or near the left branch of a v-nullcline, \(\,g^i_) \}\) for whatever inputs \((g^i_,g^i_)\) it receives, and as active when its trajectory lies on or near a right branch, \(\,g^i_) \}\). The six IN units’ projections each have subtly unique features. One point of commonality is that we have tuned the model so that each unit’s v- and h-nullclines intersect on the right branch of the v-nullcline, corresponding to some level of tonic activation, in the absence of input (cf. blue dashed v-nullclines and dotted gray h-nullclines). When receiving reciprocal inhibition from its joint core partner, however, each unit has a left branch intersection and hence cannot activate on its own, with the exception of Ret and Dep. Thus, if Pro (Lev) is active, then Ret (Dep) can nonetheless overcome the resulting inhibition and activate, suppressing its partner, but not vice versa. This tuning reflects a default bias so that the leg is on the ground, not in midair, when at rest.

The zoomed views in Fig. 6 include labels indicating the type of transition involved in the escape of each unit in the network from its silent phase. These labels can to some extent be inferred from inspection of the phase planes, but we performed a more careful, quantitative simulation experiment to solidify these determinations. Specifically, recall that when each unit first enters the silent phase, it travels along the v-nullcline \(\},0) \}\) corresponding to the inhibition that it receives from its joint core antagonistic counterpart. For each, we computed the slow timescale time of passage \(\tau\) along this nullcline from \(h=h_(0,0)\), the approximate h-value at which it enters the silent phase, to \(h=h_(},g)\), for different values of g (colored curves in Fig. 7), starting with the g value at which the saddle-node bifurcation that annihilates the unit’s silent phase critical point occurs and increasing from there. We call the computed value the clearance time and the resulting curve the clearance time curve (CTC). For our baseline rhythm, we then determined (a) the time \(\tau _\) after silent phase entry when each unit receives the excitatory input that eventually leads to its escape, (b) the strength \(g_\) of this input, and (c) the time \(\tau _\) after silent phase entry when each unit jumps up and hits the threshold \(v=-30\). We marked the points \((\tau _,g_)\) and \((\tau _,g_)\) in each plot with a plus symbol and a circle, respectively.

The locations of these points indicate the type of transition that each unit undergoes, based on the theory in Section 3.3. For FTM, the h-value for the unit must be above \(h_(},g)\) when the excitation arrives, corresponding to the plus point \((\tau _,g_)\) lying to the right of the CTC. Moreover, the difference between \(\tau _\) and \(\tau _\) must be small, corresponding to similar locations for the two marked points. This case arises uniquely only for Ext. For Pro, the first condition holds, but certain additional features are present. First, \(\tau _-\tau _\) is relatively large (as can be seen from the distance between the corresponding marked points). Second, the g value where input arrives (and where the marked points lie) is much closer to the minimal g value where the saddle-node bifurcation occurs, such that the v- and h-nullclines are likely to be much closer together. Thus, we recognize that this is a GSN transition. When \(h < h_(},g)\) at the arrival time of the excitation and hence \((\tau _,g_)\) lies to the left of the CTC, we recognize that EE has occurred. In the cases of Ret and Dep, the jump-up point is close to the CTC. For Lev and Flx, however, the distance in \(\tau\) between the CTC and the jump-up point is comparable to the distance in \(\tau _-\tau _\) for the GSN case for Pro; for Lev, the g value itself is also low and close to the saddle-node value. Thus, we denote the transitions for Lev and Flx as EE-GSN transitions, given that they have elements of EE (excitation arrives before CTC is reached) and GSN (delay between crossing the CTC and jumping to the active phase).

Fig. 7

Determination of transition mechanisms and implications for parameter robustness. Left panel: Degree of robustness of parameters associated with excitatory connections depends on the type of transition that they impact. Right panel: Determination of the transition type for the escape of each IN unit from the silent phase. For each panel, the point indicated with a plus symbol indicates the excitation level that the unit receives and the time after its silent phase entry when it arrives, while the dot is at the same excitation level but at the time after silent phase entry when it jumps to the active phase. The solid curves are the CTCs (see text)

We will now focus on the four phase transitions in the baseline rhythm: (1) activation of Ret and Flx, (2) activation of Lev and Ext, (3) activation of Pro, and (4) activation of Dep. For convenience, will refer to the marked times \(\\) from Fig. 5 to represent the times when these transitions happen. That is, we will start from the time \(\tau _1\) when Dep activates and Pro, Ext are still active, with the subsequent transition times as follows: (1) \(\tau _2\) and \(\tau _3\), (2) \(\tau _4\), (3) \(\tau _5\), (4) \(\tau _6\).

Transition (1)

When Dep activates at time \(\tau _1\), it excites both Ret and Flx. The h-value for Ret is not quite high enough to yield FTM, since \(h_ < h_(},})\), and hence the trajectory projected to the Ret phase plane undergoes an excursion on its solid green v-nullcline before Ret activates via an EE transition (e.g., at time \(\tau _2\) in Fig. 5; see Fig. 7). The h-value for Flx is approximately at the same level as the left knee of its solid green v-nullcline, \(h_ \approx h_(},})\), so any EE effect should be weak, yet the activation of Flx exhibits slightly more delay than does that of Ret (e.g., \(\tau _3\) versus \(\tau _2\) in Fig. 5). This delay arises via a GSN effect, due to the proximity of the v- and h-nullclines for Flx (Fig. 7). Thus, the activation of Flx involves a blend of the EE and GSN mechanisms, which we can also think of simply as GSN transition preceded by a brief EE period.

Transition (2)

Now, we consider the state in which Ret, Dep, and Flx are active, corresponding to the main part of the stance phase, with their joint core partners suppressed (e.g, time interval \((\tau _3,\tau _4)\), Fig. 5). In this situation, excitation signals go from Dep to Ret and Flx and from Flx to Lev and Ext. Let us consider the nullclines for the three suppressed units: Pro, Lev, and Ext. Pro is inhibited by Ret and not excited, so its green dashed v-nullcline applies and there is a stable critical point, \((v^_(0,}),h^_(0,}))\), where that intersects its h-nullcline; thus, \((E1)_e\) fails for Pro and Pro cannot activate. Although Lev is inhibited by Dep, the excitation from Flx causes its trajectory to jump from its dashed green to its solid green v-nullcline, which does not intersect its h-nullcline. Thus, \((E1)_e\) holds for Lev and this represents an EE configuration (Section 3.3). After some time of transit up the left branch of this v-nullcline, at time \(\tau _4\), Lev reaches \((v^_(},}),h^_(},}))\), experiences an additional GSN effect, and then fully activates, suppressing Dep, yielding another EE-GSN transition overall (Fig. 7).

Ext also receives excitation, along with inhibition, from Flx, so why does Ext not lead this transition? The solid green nullcline for Ext in Fig. 6 corresponds to excitation from Lev to Ext, not from Flx to Ext, which is much weaker. The excitation from Flx pushes the Ext trajectory just slightly below its purely inhibited (dashed green) v-nullcline. Although it cannot be definitively seen in Fig. 6, Ext still has a stable critical point, \((v^_(},}),h^_(},}))\), for this level of excitation, violating \((E1)_e\). When Lev activates, it excites Ext, the Ext v-nullcline drops down to the solid green location, with \(h_> h^_(}+},})\), and Ext immediately follows Lev into the active phase via FTM (Somers and Kopell, 1993; Rubin and Terman, 2002) (Section 3.3; Fig. 7).

Transition (3)

Once Lev and Ext activate, Lev excites Pro. So why does Pro not activate together with these other units? Notice in Fig. 6 that the projection of the trajectory to the Pro phase plane does jump away from the dashed green v-nullcline when this excitation arrives, at \(t=\tau _4\). Its solid green v-nullcline, corresponding to its being excited (by Lev) and inhibited (by Ret) does not intersect its h-nullcline; that is, the onset of excitation to Pro from Lev induces a saddle-node bifurcation of equilibria for Pro and \((E1)_e\) holds. Nonetheless, the two nullclines lie extremely close together. Their proximity yields a region in phase space where both v and h change slowly for the Pro unit, with \(h_ \approx h_^(},})\). The resulting GSN (Section 3.3; Fig. 7) delays the onset of Pro activation, such as from time \(\tau _4\) to time \(\tau _5\) in Fig. 5.

Transition (4)

Just as the activation of Lev excites Pro, the activation of Ext excites Dep. Dep does not have a GSN and \((E1)_e\) holds, so why does it not become active with Ext? Recall that Dep first left the active phase in Transition (2), and the value of its h-coordinate is low when Ext activates (e.g., at time \(\tau _4\) in Fig. 5). Thus, even though the projection of the trajectory to the Dep phase plane lies on the solid, excited v-nullcline of the Dep unit’s green (inhibited) nullcline pair, which does not intersect its h-nullcline, it has \(h_ \ll h_^(},})\) and the projected trajectory must undergo a long excursion up the nullcline’s left branch before reaching its knee. When this excursion is completed, Dep activates (e.g., at times \(\tau _1\) and \(\tau _6\) in Fig. 5), representing another EE transition (Fig. 7).

In summary, we see that the transitions in which units activate during this multi-phase baseline rhythm feature all three of the mechanisms that we presented previously (Section 3.3). Of these, both EE and GSN mechanisms lead to a delay between the activation of one unit and the activation of another unit that it excites, while FTM does not.

The three mechanisms yield distinct predictions about sensitivity of activation timing to parameter variations (see also Section 3.3). In the GSN mechanism, the delay time would be expected to be sensitive to small changes in the v-nullcline position, which could arise due to changes in excitation strength or in certain parameter values intrinsic to the unit’s model equations. Parameter changes that push the v-nullcline and specifically \(h_\) even in a little bit lower would be expected to significantly shorten the delay in activation and could switch ghost transitions to FTM.

In EE, we expect more robustness to parameter variations as long as the associated trajectory’s jump to the excited v-nullcline’s left branch happens far from its left knee, such that the slow excursion on the left branch is preserved as parameters are varied. For excitation onsets relatively later within a unit’s silent phase, however, a trajectory’s h-value will be closer to that of the relevant left knee, such that parameter changes that affect nullcline position may switch a transition from EE, with \(h < h_\) and some activation delay, to FTM, with \(h> h_\) and a much shorter delay.

Finally, FTM itself occurs abruptly on the fast timescale, so once a transition is within the FTM regime, activation timing should be robust to subsequent parameter changes unless they are so extreme that they switch the transition mechanism from FTM to one of the other types. We will explore issues of robustness in much more detail in Section 4.3.

4.3 Changes to different parameters have different effects on rhythms

To find parameter sets that would yield the baseline rhythm, we tuned the three parameter classes \(\\}\), \(\}\}\), and \(\}\}\) individually across units and coupled unit pairs. Here, we consider the robustness of the rhythm to various forms of changes in the values of these key parameters. We note that we are considering a model with a highly reduced representation of the signaling associated with feedback. Potentially, more sophisticated modeling of these signals and interactions across limbs could promote robustness in a more biological setting; nonetheless, we assume that maintenance of a rhythm satisfying (P) across some range of parameter variations represents an important feature for our model. We will consider four issues in the context of robustness: (1) How robust is the baseline rhythm to changes in different classes of parameters, and why does robustness vary across these classes and specific elements within them? (2) Which parameters can be tuned to adjust rhythm frequency, and which rhythm features vary under this tuning? (3) When (P) is sensitive to changes in certain parameters, to what extent can variation of other parameters compensate to maintain (P) (cf. Marder et al. 2015)? (4) Can inhibitory interactions across joint segments enhance robustness?

4.3.1 Parameter sensitivity relates to the transition mechanisms affected

First, we tested the effect of multiplying all values in a selected parameter class by \((1+\epsilon )\) as we varied the scaling parameter \(\epsilon\) over a small interval around 0. We found that the most sensitive parameter class was \(\\}\), which could only tolerate scaling by \(\mathcal (10^)\) values of \(|\epsilon |\).

Although this uniform scaling was quick to implement as a preliminary step, and could represent biological variability in some factor that is shared across neurons or synapses, we would expect that variability in general could occur at the level of individual neurons and synapses. Thus, we next turned to variation of each individual element of each of the three parameter classes. In our simulations, we tested for robustness with respect to variations of up to \(\pm 5\%\) in parameter values. We expected to find the least robustness with respect to parameters that affect the activation of Pro, since Pro activation relies most heavily on the ghost transition mechanism, which is sensitive to even small changes in v-nullcline positions. Flx and Lev activations also feature a weaker form of the GSN effect, and hence seemed like another natural source of parameter sensitivity.

Overall, we found the most robustness with respect to variations in the excitatory, inter-segment \(\}\}\) parameter class. The network continued to produce a rhythm that satisfied (P) over \(\pm 5\%\) variation of all such parameters except \(}\), where robustness was limited, and \(}\), which failed (P) specifically for \(+5\%\) variation. Since \(}\) affects the position of the solid green Pro v-nullcline that is critical to Pro activation (Fig. 6), non-robustness with respect to \(}\) matches our expectations. The much milder non-robustness with respect to \(}\) is also consistent with our expectations based on the GSN mechanism of Lev.

Robustness with respect to changes in \(\\}\) and the inhibitory, intra-joint coupling \(\}\}\) was much more limited. Within each of these classes, however, certain parameters yielded more robustness than others. The parameter \(g_\) was the most sensitive of the external drive parameters, tolerating variations of less than \(\pm 1\%\), again consistent with the sensitivity of ghost transitions. As this parameter was varied, the most fragile property of (P) was the synchronized inactivation of Pro and Ext. Decreases in \(g_\) delayed the activation of Pro and the inactivation of Ret. This change propagated into delays in subsequent Ret activation, which caused Pro to fall silent too long after Ext, and worsened across successive cycles. Similarly, increases in \(g_\) had the opposite effect, with successively earlier Pro activations and Ret activations, the latter of which desynchronized the inactivation of Pro and Ext.

We also found quite limited robustness in \(g_\) and, to a slightly lesser extent, in \(g_\). Changes in \(g_\) led to similar issues with desychronization of Pro and Ext inactivation, which we consistently found to be a vulnerable aspect of the rhythm, stemming from small changes in the timing of Pro activation through the GSN mechanism together with the lack of any direct connection to help coordinate Pro and Ext activity (see Fig. 1). Decreases of \(g_\) also harmed the relative timing of Pro and Ext inactivation: with small decreases, Pro inactivated too early relative to Ext, and with slightly larger drops it also inactivated too soon after Lev inactivation. Lower \(g_\) delays Lev’s escape via its GSN, which correspondingly delays its recruitment of Ext by FTM, but Pro and Ret are less affected, which throws off the relative timing in the rhythm. On the other hand, increases in \(g_\) caused the rhythm to fall apart entirely. These increases allowed for earlier Lev escape. As a result, the excitation from Lev to Pro arrived when \(h_\) was smaller, making it more difficult for Pro to undergo its GSN escape. This problem worsened across cycles, until eventually Pro failed to activate on one or more cycles.

The components of \(\}\}\) are mostly sensitive, with the exception of \(}\) and \(}\); We expect that the stronger sensitivity to \(}\) than \(}\) comes from the fact that the excitation level specifically affects the activation process. Although certain units are excited throughout their active phases, voltage values lie near the excitatory reversal potential during these phases, so small changes in excitation levels do not have much effect there. The inhibition to a unit from its joint partner, however, affects its evolution throughout its entire silent phase, potentially leading to stronger impact.

4.3.2 Parameter sensitivity can be mitigated by variation of compensating parameters

Although heterogeneity across individuals, neurons, and synapses could naturally occur at the individual parameter level, experiments have also revealed evidence of co-variation between values of different parameters (Marder, 2011; Lamb and Calabrese, 2013; Goaillard and Marder, 2021). Moreover, there are biological mechanisms for which variations would translate into changes in multiple parameter values; for example, variation of the density of the synaptic neurotransmitter receptors on a post-synaptic neuron would affect the strengths of all inputs to that neuron that use that neurotransmitter, while variation of the rate of synaptic release by a pre-synaptic neuron would affect the levels of its coupling to all of its synaptic targets. Such joint changes can occur biologically through the action of neuromodulators (Marder and Bucher, 2001; Svensson et al., 2001; MacKay-Lyons, 2002). Therefore, we next considered the extent to which coordinated changes in the values of two or more parameters in the model could improve the robustness of its output pattern.

To start with, we considered the sensitive parameter \(}\). The model maintained a rhythm that satisfied (P) with nearly \(50\%\) cuts in \(}\) when these were accompanied by small, proportional increases in \(g_\). These compensations maintained the ghost mechanism as needed for (P). We could also maintain a (P) rhythm with up to \(10\%\) increases in \(}\) by matching each increase with a larger, proportional decrease in \(g_\). With stronger \(}\), Pro activates earlier and inactivates inappropriately before Ext. A weaker \(g_\) delays escape by Ret and hence leads to a prolonged period of Pro activation, which maintains the desired timing of Pro relative to Ext. Despite the shorter delay from Lev activation to Pro activation in this regime, these changes preserve (P), as long as they are not too large. See Fig. 14 in Appendix B for a numerical example.

In a similar vein, we found that we could achieve robustness with respect to at least \(\pm 5\%\) variation of every one of the \(g_\) parameters through compensating variation of another 1-2 parameters. Variations in all the \(g_\) parameters, except for \(g_\), could be fully compensated by either introducing stronger, same-sign variations in \(}\) or by opposite changes in interjoint excitatory coupling \(}\). For instance, variations in \(g_\) could be fully compensated by corresponding, opposite changes in \(}\), based on similar reasoning relating to the Pro active period, or by introducing stronger, same-sign variations in \(}\). Among the most sensitive parameters, compensation for \(g_\) is achieved by increasing both \(}\) and \(}\), while introducing stronger, same-sign variations in \(}\) together with same-sign variations in \(}\) yielded robustness to \(\pm 5\%\) variations in \(g_\).

Analogously, while almost all \(}\) parameters are sensitive, they can tolerate greater variations if they are varied in coordination with \(g_\) parameters. These results highlight the intrinsic duality between \(}\) and \(g_\): changes in one can be systematically offset by adjustments in the other, reinforcing their compensatory relationship.

4.3.3 Neuromodulation can tune rhythm frequency and relative phase durations

To complement this analysis, we specifically considered sets of parameter changes corresponding to groups of synapses with common sources or common targets, which could most naturally covary – or be co-varied – through biological mechanisms such as neuromodulation. For example, jointly scaling \(}\) and \(}\) by a common factor would correspond to modulating Dep outputs, and this manipulation maintained a rhythm satisfying (P). We found a robustness to combined variations of strengths of excitatory outputs from most individual sources and to combined variations of strengths of all excitatory inputs to most shared targets. Certain exceptions occurred that were, not surprisingly, associated with the Pro unit. The parameter \(}\) is sensitive, and in order to make rhythms robust to joint scaling of \(}\) and \(}\), the two excitatory outputs from Lev, we had to adjust either \(g_\) or \(g_\) as well. As another example, the only excitatory input to Pro is \(}\), which we have already seen is highly sensitive to variations.

One interesting case that we highlight is joint scaling of all excitatory outputs from Flx, namely scaling \(}, }, }\) by a common factor \(\lambda\). We could achieve more than a 10\(\%\) variation in the period of a rhythm that satisfies (P) using scaling factors of around \(\pm 5\%\), i.e., \(\lambda \in [0.95, 1.05]\). These changes in period arose through changes in the activation periods of Ret, Dep, and Flx (see Fig. 8), associated with the stance phase of the locomotor pattern, which is the aspect of the locomotor rhythm that varies under biologically observed period changes as well (Gabriel and Büschges, 2007; von Uckermann and Büschges, 2009). See Fig. 15 in the appendix for time series illustrations of Pro and Ret, along with changes in period, swing, and stance durations at extreme values of \(\lambda\) (i.e., \(\lambda = 0.95\) and \(\lambda = 1.05\)).

Upon further investigation, we found that joint scaling of all the inputs to Lev, namely scaling \(g_\), \(}\), \(}\) by a common factor \(\lambda \in [0.95, 1.05]\), achieves the same strong robustness of (P) and variation in period. In fact, the key connection parameter underlying this effect as well as that from scaling excitatory Flx outputs is \(}\). Note that when Flx activates, it excites Lev and Ext. The next unit to escape after Flx is Lev, and this escape occurs through the EE mechanism associated with the excitation from Flx, with the actual moment of escape happening when the trajectory reaches \(h_ = h_(},})\). The h-coordinate of this knee is modulated by \(}\); stronger excitation leads to a lower \(h_\) and hence earlier escape, which corresponds to a shorter silent phase of Lev and active phase of Dep (Fig. 8, stance panel) and hence contributes to the shortening of period. Once Lev escapes, it excites Ext, which escapes via FTM, a highly robust mechanism, leading to a shorter active phase of Flx, and it also excites Pro, which shortens the Ret active phase, both of which further decrease the period. On the other hand, the fact that Lev, Ext and Pro jump up with different timing does not affect the h values where their respective joint block partners, Dep, Flx, and Ret, return to the silent phase, since in all cases, that happens very close to a critical point on the right branch of the corresponding v-nullcline, associated with condition (E3). Thus, Dep, Flx, and Ret will not undergo changes in their silent phase durations before their next escape, and correspondingly there are no changes in the Lev, Ext, and Pro active phase durations (Fig. 8, swing panel). In fact, there is no other parameter that we have identified that can induce such a robust variation of cycle period, because of the unique role of the excitation from Flx to Lev. Indeed, although Dep and Ret both activate via an early excitation mechanism i