# Toward robust phase-locking in Melibe swim central pattern generator models

## Abstract

Small groups of interneurons, abbreviated by CPG for central pattern generators, are arranged into neural networks to generate a variety of core bursting rhythms with specific phase-locked states, on distinct time scales, that govern vital motor behaviors in invertebrates such as chewing, swimming, etc. These movements in lower level animals mimic motions of organs in higher animals due to evolutionarily conserved mechanisms. Hence, various neurological diseases can be linked to abnormal movement of body parts that are regulated by a malfunctioning CPG. In this paper, we, being inspired by recent experimental studies of neuronal activity patterns recorded from a swimming motion CPG of the sea slug Melibe leonina, examine a mathematical model of a 4-cell network that can plausibly and stably underlie the observed bursting rhythm. We develop a dynamical systems framework for explaining the existence and robustness of phase-locked states in activity patterns produced by the modeled CPGs. The proposed tools can be used for identifying core components for other CPG networks with reliable bursting outcomes and specific phase relationships between the interneurons. Our findings can be employed for identifying or implementing the conditions for normal and pathological functioning of basic CPGs of animals and artificially intelligent prosthetics that can regulate various movements.

Many abnormal neurological phenomena are perturbations of normal functions of the underlying mechanisms governing the animal behaviors, specifically movements. Repetitive behaviors are often associated hypothetically with the phenomenon of rhythmogenesis in small networks that are able autonomously to generate or continue, after induction a variety of activity patterns without further external input, abrupt or not. The goal of this modeling study is to identify decisive components of a biologically based CPG that has been linked to a specific motion in a lower order animal, sea slug Melibe leonina, which produces a specific bursting pattern, as well as to identify components that ensure the pattern’s robustness.
Due to the recurrent nature of such bursting patterns of self-sustained activity we employ Poincare return maps defined on phases and phase-lags between burst initiations in the interneurons to study quantitative and qualitative properties of CPG rhythms and corresponding attractors.

The proposed approach is specifically tailored for various studies of neural networks in neuroscience, computational and experimental. Development of such tools and our understanding of such CPGs can be applied to gain insight into governing principles of neurological phenomena in higher order animals and can aid in treating anomalies associated with neurological disorders.

## I Introduction

A central pattern generator (CPG) is a neural network of a small group of synaptically coupled interneurons that is able to generate single or multiple rhythmic outcomes without external [sensory] feedback (1); (2); (3). CPGs establish and govern various motor behaviors of animals such as swimming, crawling, walking etc (4); (5). In addition, the mechanisms of such motions are evolutionarily conserved and can be related to rhythmic motions of various body parts, such as heart, lungs, legs etc, of higher order animals. Such rhythmic outcomes, often viewed as bursting patterns, can be differentiated by several timing properties, including specific and robust phase-locked states between well orchestrated interneurons within the specific CPG of a particular animal. As such, the behavior controlled by the CPG can be disrupted, or halted, after a component neuron of the network is blocked or temporarily inhibited. This would indicate that the rhythmic outcome results from synergetic interactions of all contributing interneurons, which may not be necessarily endogenous bursters, when isolated from the network. The robustness of a rhythmic outcome is an essential property allowing the CPG to withstand or recover from perturbations, a lack of which could be the expression of various neurological diseases and disorders.

Identification and modeling of a CPG underlying an animal behavior is a real challenge due to a number of factors. The realization of a behavior may require components, other than interneurons, such as synapses, which can be fast and slow, inhibitory, excitatory, or electrical, etc. Furthermore, interneurons, which are networked within the CPG for one behavior, may contribute to another behavior as well, i.e. be multifunctional (6); (7). A whole CPG can be also multifunctional if it governs more than one behavior, in contrast to a dedicated CPG which is arranged for the purpose of a single locomotion. Modeling studies, mathematical and computational, have proven to be useful for gaining insights into operational principles of CPGs (8); (9); (10); (11); (12); (13), in particular, multifunctional ones (14); (15); (16). This study has been inspired by recent experimental studies of neuronal activity patterns recorded from the identified interneurons of a [dedicated] CPG governing the swim pattern of a sea slug Melibe leonina (17); (18); (19) and possible consequences of understanding neurological phenomena in general.

The occurrence of phase-locked bursting patterns is naturally observed in the voltage traces simultaneously recorded from a few interneurons, which are argued to belong to the same CPG. In such recordings, specific delays between initiations of the active (spiking) phases of the bursting interneurons are well maintained. It is unclear so far whether the delays or phase-lags between bursting interneurons are essential for the CPG. An argument supporting the latter hypothesis is that the same phase-lags between burst initiations have been recorded in both juvenile and adult animals (18), while the frequency and hence delays between burst initiations in voltage traces can vary substantially during animals’ life spans. Bursting neurons display multiple action potentials during their active phases of bursting and remain hyper-polarized during their inactive phases. The bursts of action potentials or spikes correlate with neurotransmitters’ release that allows the neurons to interact. Hence, the specific delays between the bursting patterns can be meaningful to and explanatory for CPG formation mechanisms (20); (21); (22); (23); (24); (25); (26).

The question of whether a neuron belongs to a CPG is not easy. Potentially it can be resolved by linking the sequence of bursting delays with some timescale of movements during the behavior. An alternative is a perturbation approach when a targeted neuron is temporarily injected with a polarized current to trace down a metamorphic effect or a qualitative change on the network. A modeling study of the CPG, however, remains a practical approach for singling out particular features among networks that are vital for its own proper functioning (35); (36). In this paper, we will start with the consideration of rather idealized, symmetric configurations of CPG networks whose repetitive bursting patterns or rhythmic outcomes are known a priori. By doing so, we will probe tools developed for uncharacteristic dynamics, including a possible co-existence of several patterns in other CPG configurations. Next, we introduce, gradually, variations in the parameters to match our findings with the recorded electrophysiology of the animal. The strategy should allow us to find a minimal wiring for synaptic connections that gives rise to (robust) neuronal dynamics observed in experimental studies.

We begin by introducing the dynamical toolkit in the Methods section that is used throughout this study. It is followed by a section on a CPG example made of uncoupled half-center oscillators (HCOs). A HCO is composed of two bilaterally symmetric cells reciprocally inhibiting each other to produce alternating (anti-phase) bursting patterns. Such a pair can burst in-phase too (33); (34), when, for example, it is exogenously driven by a pre-synaptic interneuron of the network (12). We will consider two cases of uncoupled cs, homogeneous and heterogeneous. By comparing their dynamics, we identify the conditions giving rise to robust and unique phase-locked pattern formations, such as ones recorded in the experiments. As a next step, we introduce additional synaptic connections, which are arguably known to exist in the circuitry model of the swim CPG and study their roles in regulating the mathematical models of CPG. Then, we can evaluate the parameter range for the network heterogeneity, which sustains the plausible phase-locked patterns, and show how the latter depend on the network configurations. Finally, we construct return maps for the phase-lags between both HCOs, not interneurons. This will further reduce the original problem (coupled 12 ordinary differential equations) to low-order maps to study synchrony, stability, coexistence and bifurcations of the bursting patterns. In the Discussion section we will also address the challenges and future directions and application for the proposed dynamical framework.

## Ii Methods

In this study of 4-cell CPG networks, we employ a generic Hodgkin-Huxley-like model of an endogenously bursting interneuron as an elementary block; the building blocks will be the HCO that can intrinsically burst either in-phase, in general, or anti-phase, in particular. This reduced model, describing the dynamics of a leech heart interneuron, has been extensively studied and biophysically calibrated to demonstrate a variety of activity patterns typical for various invertebrates (28); (29). Depending on external drive, the interneuron model can produce tonic spiking activity, be an endogenous burster, or settle down to hyper- and depolarized quiescence states. The model has turned out to be very reach dynamically as can demonstrate a number of global bifurcations at the activity transitions, like the blue sky catastrophe, bi-stability and chaos due to homoclinic saddle orbits (30). In this study, individually each post-synaptic interneuron is an endogenous burster that can become temporarily shut at the hyper-polarized state by an inhibitory current originating from pre-synaptic interneuron(s) (12); (15).

This level of accuracy is important for understanding CPG mechanisms as models must be compared with real animal behaviors for testing our hypotheses, however, eventually model parameters would be modified to fit neurological phenomena in mammals, in particular humans, for investigating disorders with neurological origins.

Below, we will consider several types of CPGs made of the interneuron models weakly coupled by synapses, chemical: inhibitory and excitatory, and electrical, referred to as a gap junction. The equations of the coupled model are given in Appendix. Chemical synapses, inhibitory and excitatory, are described within the framework of the fast threshold modulation (FTM) paradigm (32), which has been proven to meet some basic conditions for coupled bursters (33); (34). The strength of coupling is controlled by the maximal conductance, , for the synaptic current. Besides , as its magnitude should be sufficient to guarantee a slow rate of progression of bursting patterns, transitioning toward a phase-locked state, if any. We ensure that the convergence is not due to a symmetry of network interactions; some deviations, , from the nominal values are introduced in the inhibitory synapses: . Unless otherwise mentioned, , , , , and .

We must point out that a transient trajectory can converge to a network attractor rather quickly even in a weakly coupled case (). Such a quick convergence can occur when the endogenously bursting interneuron is initially close to a transition to a quiescent steady state through a slow-time scale bifurcation like saddle-node or homoclinic. In this case, the trajectory can come close by or cross back and forth the corresponding (bifurcation) boundary when an interneuron receives (or is released from) a flux of inhibition from another pre-synaptic interneuron on the network. Physiologically speaking, such neuromodulation can be viewed as an analogues to the mathematical phenomena of bifurcations through perturbations. Natural substances such as serotonin released by the animal can alter intrinsic properties of the individual neurons affecting the efficiency of a CPG and vary its temporal characteristics without breaking the bursting pattern per se. (31).

In the CPG mathematical model we set all the parameters so that the individual and networked interneurons remain endogenous bursters. The duty cycle of bursting of the interneurons, which is a fraction of the period during which the interneuron is active, persistently stays around , i.e. the burst (active) durations and hyper-polarized (inactive) periods are almost equal. Figure 1 shows a typical bursting pattern that resembles the experimental recordings from four interneurons of the Melibe swim CPG.

The specific delays between the burst initiations between the interneurons are the key characteristics of each given CPG. The core idea underlying our computational tools is inspired by “wet lab” experimental observations and therefore tailored for neuroscience (15). In essence, it requires only the voltage recordings from the mathematical interneurons, and therefore does not explicitly rely on the gating variables from a Hodgkin-Huxley type model. We intentionally choose the phases based on the membrane voltages, as these are basically the only variables that can be experimentally assessed and measured. Moreover, as in wet experiments, we have control over, and hence can maintain the initial delays, or phase distribution by releasing the interneurons from inhibition at various times after or prior of the release of the reference neuron.

The phase relationships between the coupled interneurons are defined through specific events, , that occur when their voltages reach an auxiliary threshold, V, set above the hyperpolarized voltage and below the spike oscillations. Such events indicate the initiation of the sequential bursts in the interneurons, see Fig. 1.

We define a sequence of phase-lags through the delays in burst initiations relative to that of the reference neuron 1, normalized over the current network period, or, specifically, the burst recurrent times for the reference interneuron, as follows:

(1) |

An ordered triple, , defines a forward iterate, or a phase point (see Fig. 2), of the Poincaré return map for the phase-lags: . A sequence, , yields a forward
*phase-lag trajectory*, , of the Poincaré return map on a 3D torus with phases defined on modulo 1, see Fig. 4.

Based on the experimental recordings of interneurons of the Melibe swim CPG the authors (18) suggest a possible architecture for its network. Its wiring schematics shown in 2 includes two core half-center oscillators as the network building blocks: HCO1 (top, shown in blue) and HCO2 (bottom, shown in pink). The interneurons of each HCO burst robustly in anti-phase while the animal swims. The interneurons within a HCO are known to inhibit contralaterally each other: — and —, while 1 and 2 excite — 3 and 4 ipsilaterally, resp. In addition there is an electrical coupling through a gap junction between the top interneurons, 1 and 2. Strong electrical ipsilateral coupling between 1 and 3, as well as 2 and 4 on the other side, makes them oscillate together.

Traces of the bursting membrane potentials of the interneurons of the CPG are used to derive the phase-lags according to the method illustrated in Fig. 1. By varying the delays between burst initiation in the reference interneurons and release times of three other ones from inhibition, we obtain a dense array of initial phase distributions. Then, the tuples of three phase-lags are recorded at every cycle of the network bursting as the time progresses. Figure 2 illustrates a typical evolution of the resulting sequences, , and plotted against the burst cycle number, , converging to a phase-locked state around ; (unless otherwise mentioned, green, black and blue are the color-codes for the phase-lags, respectively.)

We represent the sequence as a forward trajectory on a 3D-torus as the burst cycle number, , is increased. Geometrically, the 3D is viewed as a solid unit cube shown in Fig. 2. The opposite sides must be identified due to the cyclic nature of the phase; this implies that trajectory leaving the cube through one of its six sides, will be wrapped around to re-enter through its opposite side and so forth.

Alternatively, we can consider progressions of the phase lags individually, in terms of 1D Poincaré return maps: (Fig.2). Due to weak coupling, sequential iterates of the maps do not jump far apart from each other. This allows for having them connected into “continuous” trajectories, in order to follow their evolution in forward time. The green sphere, in Fig. 2, representing the initial phase-lags tuple, corresponds to the beginning of a trajectory. Such initial phase-lags are uniformly distributed on a lattice within a unit cube, see Fig. 2).

The tuple, , of the phase-lags represents the state of the network at the -th burst cycle, because it captures the temporal bursting activity of all four neurons of the CPG. In what follows, we will explore visually how the network state progresses by following the evolutions of various initial phase-lags, which can converge to a single or multiple attractors. Such an attractor of the map corresponds to a stable bursting pattern configurations of the CPGs in question.

The initial distribution of network states or phase-lag tuples, sampled uniformly to form the cubic lattice within the cube, will shift toward attractors as the burst number increased. In the case when such an attractor is a fixed point, its coordinate corresponds to stable phase-locked state of the bursting pattern with specific time delays between the interneurons of the CPG. It is equivalent that a single stable fixed point of the return map describes a single robust pattern of the dedicated CPG, as all initial delays will ultimately lead to the same bursting pattern. As the parameters of the CPG are changed, the stable fixed point can bifurcate, for example vanish or loose the stability, thus giving rise to another attractor such as an invariant circle. In the latter case, the phase of the bursting pattern are no longer locked, but vary periodically, or even show some aperiodic, chaotic dynamics.

If the corresponding return map shows two (or more) attractors, then depending on initial phases the multifunctional CPG can respectively produce several patterns. Perturbations, such as noise or external polarized currents, can cause sudden or unforeseen jumps between the attractors, resulting in switching between the corresponding patterns, for example between in-phase and out-of-phase bursting.

It is worth noticing that in weakly coupled cases, one should consider constructing 1D Poincaré maps for every -th phase-lags. This procedure gives a map of the degree , with a “flatter,” so to speak, graph adjoining the -line at a stable fixed point. Figure 2 displays such maps of degrees 1, 5 and 10, with the corresponding fixed points around , and for the individual phase-lags as predicted from trace progressions depicted in Fig. 2. Note that the fixed points are approached from one side only.

It is our working hypothesis that a stable state of the CPG is defined by interplays of synaptic strengths, rather than by specific wiring of the network. While wiring can be a necessary condition for the network to produce bursting pattern(s) , its configuration does not provide the sufficient condition for the robustness of the latter. Thus, the problem of the robustness of the pattern is reduced to the stability conditions and identification of bifurcations of the corresponding fixed points of the map for the phase-lags. The current state is such that the maps have to be visualized to identify and classify all patterns resulting from different initial phase relations between the four bursting interneurons.

The following concern must be also addressed: how can changes in phase-lags relate to the periods of the two HCOs forming the 4-neuron CPG. Namely, whether it is crucial that the periods of both individual HCOs are not the same, and how coupling affects the period of the whole network. Figure 3 addresses this issue: let HCO2 have the period, , slightly longer (shorter) than the period, of HCO1 (i.e. or ), then phase-lags between them, here , will increase (decrease) in the CPG, which begins from the same initial conditions. This relation between the phase-lags and the periods is discussed in Appendix.

## Iii Results

The main question that we aim to address in our modeling study is: what synaptic connections are the key ones that lead to the experimentally observed (i.e. stable) phase-locked bursting in the CPG.

Let us first consider a network configuration of the CPG with only contralateral inhibitory connections of the strength , which is a half of that of the strongest inhibitory synapses in the CPG (). Later, we introduce the ipsilateral excitatory synapses ( and ), followed by the electrical synapses or gap junction () between the interneurons 1 and 2 to examine transformations of bursting patterns through bifurcations, if any, of attractors in the corresponding 3D maps on the torus.

We have performed comprehensive simulations and further visualization of solutions of return maps for the sequences of phase-lags for the CPG configurations. For the given CPG model, the corresponding 3D Poincaré map is shown in Fig. 4. It displays multiple transients converging to a stable fixed point corresponding to the phase-locked lags within the bursting pattern. In this figure, a green dot indicates the beginning point of a transient of the map in every simulation run. By releasing trajectories from a dense, homogeneously distributed grid of initial phases (conditions) spread over the bursting periodic orbits, one can obtain a complete portrait of the phase space of the return map on the unit cube. This 3D portrait in Fig. 4 shows attractors, separating saddles and invariant subspaces. Here, the red circles indicate the locations of saddles, or turning points, while the cyan and the blue circles correspond to the steady states of the network, i.e. stable fixed points of the map.

The fixed point represented by the solid cyan circle in Fig. 4 has the following coordinates: and . This means that the interneuron 1 and 3 burst in-phase with respect to each other, while they keep bursting in anti-phase with their counterparts: interneurons 2 and 4. In other words, the CPG is made of the HCOs bursting in-phase, synchronously, with respect to each other. By inspecting the stability of the fixed point in restriction to the side, of the cube, one can conclude that there are two ways which give rise to this bursting rhythm. The independent HCO2 can either slow down, or catch up, to line up with the inhibitorily driven HCO1. One can conclude here that, in essence, the phase relations between the interneurons of this CPG is effectively reduced to that between the HCOs, provided that both maintain anti-phase bursting endogenously.

There is another stable fixed point of a smaller basin in the phase space of the map for this CPG configuration. Its location is indicated by the solid blue circle, , of the side, given by , of the cube. The coordinates of this fixed point correspond to a rather constrained rhythm: the in-phase bursting interneurons 1 and 2 of HCO1 are driven, causing a phase-lag, by the also in-phase bursting interneurons 3 and 4 of HCO2. Therefore, near this fixed point the 4-cell CPG acts as a 2-cell network, which is equivalent to one interneuron inhibiting the other with a double drive.

Some projections of the map on the unit cube can be misleading for evaluations of the locations of the stable fixed points, as one has to see the orthogonal projections, along with the phase-lag progressions plotted against the bursting cycle. Note that interpretations of such progressions (projections too) could be a challenge in a multi-stable case where overlapping trajectories tend to converge to several attractors. To give a comprehensive overview of the dynamics of the 3D map we will also utilize such orthogonal 2D projections and frequency coint distributions to understand better the behavior of its transients. Of particular interest are how they converge to the attractors, as in some cases the convergence can be achieved only from one side. In terms of the map for a dedicated CPG, this means strengthening the stability of the single fixed point from any direction in the unit cube, which becomes its global attraction basin. Addition of synapses to the CPG may make it multifunctional, so it is imperative to know how this changes quantitatively the number of attractors and qualitatively the stability conditions for the bursting patterns.

We start the next section with a discussion on a dissected CPG with uncoupled HCOs. While this may seem trivial, it should give us a reference framework necessary for singling out the underlying organizations of the phase-lag trajectories resulting from addition of basic synaptic connections. In the 4-cell network, the connections, contralateral inhibitory and ipsilateral excitatory, should promote, not conflict, the robustness of bursting outcomes of the CPG.

### iii.1 CPG dissection in homo- and heterogeneous HCOs

A network state will correspond to a behavior pattern if it evolves on par with the behavior. By introducing variations in carefully chosen synaptic connections of the network, we make predictions and match outcomes with the expected behaviors, in order to identify the CPG mechanism. In this study, we assume that a persistent phase-locked state underlies the Melibe swimming behavior. While the ideal mathematical model must reflect all experimentally observed features of the biological CPG, a reduced model is intended to describe only some likely mechanisms giving rise to stable phase-locked bursting patterns such as the one depicted in Fig. 1.

In order to elucidate how CPG networks operate in general, and, in particular, how the Melibe CPG, robustly produces the single pattern with the constant phase-lag, we apply a bottom-up approach. This approach is used for identifying and differentiating the features that persist as the network configuration becomes more plausible in comparison with the biological CPG architecture. For example, one pair of uncoupled HCOs suffices to produce anti-phase bursting patterns observed in the voltage traces recorded from four interneurons. In addition, the capacity of pattern generation of 3-neuron motifs made of reciprocally inhibitory interneurons is well understood (7); (15).

It was shown recently that under certain conditions, fast non-delayed reciprocal inhibition within a stand-alone pair of similar neurons may lead to synchronous, in-phase bursting (33); (34). So, for the sake of generality, we set the parameters of the individual interneurons and the cross-coupling some different to guarantee that anti-phase bursting is the only stable pattern in either HCO (12); (15).

Figure 5 recaps some findings for the dissected CPG made of uncoupled () and homogeneous HCOs: all maximal conductances are equal . Figure 5 shows a few samples of the phase-lags progression, , and , plotted against the cycle number. It shows that the network transients converge to more than a single attractor. Observe too that convergence rates to the phase-locked states are predictably equal in this homogeneous case. Next to it is the frequency count distribution of 448 initial network states after 100 burst cycles. The diagram depicts two dominating peaks in green for () corresponding, respectively, to in-phase and anti-phase bursting in HCO1 (HCO2). As the HCOs are uncoupled, the distributions for (in black) and (in blue) are uniform. Since the HCOs receive no inputs from each other, they evolve independently resulting in the phase-lags between the reference neurons in each to be arbitrary.

Figure 5 represents the 2D -projection of the phase-lag transients. In it, green dots unmask a uniformly distributed lattice of initial states of the network, while blue crosses mark the ends of the phase-lag transients. Note again that when, a phase-lag goes over 1, its value is reset by modulo 1. Black squares in Fig. 5 indicate regions of rather slow evolution of the network transients, where sum of the phase-lags per cycle shift less than over the last ten cycles.

The 3D phase space of the unit cube for the uncoupled HCOs is given in Fig. 5. It shows two nodal attractors corresponding to the anti-phase bursting in each HCO. Because of the equal coupling weights, each node is “symmetric” thus indicating even convergence rates in both HCOs. The lack of interaction between HCO explains the presence of an invariant curve, given by the constraints (mod 1) and . This invariant line corresponds to anti-phase bursting in the uncoupled HCOs. The curve is stable in directions transverse to it, and neutrally stable along it, i.e. the phase-lags on it do not shift in the homogeneous case. In other words, each HCO tends typically, for most initial conditions, to anti-phase bursting.

Next, let us consider a heterogeneous CPG such that the reciprocal inhibitions in HCO1 are two times less than those in HCO2: , where . The corresponding return mappings are shown in Fig. 6. Symmetric and phase-lag projection shows the persistent nodal attractor corresponding to the endogenously anti-phase bursting HCOs. The quantitative changes of the heterogeneous case compared to the outcomes of homogeneous one is that the fixed point has the leading (horizontal) and strongly stable (vertical) directions (Fig. 6) due to, correspondingly, faster and slow convergence rates to the anti-phase bursting in HCO2 and HCO1 with stronger and weaker reciprocally inhibitory synapses. As a result, transients, which have converged to the invariant line, (mod 1) (Fig. 6) slide slowly along it. The slow speed rate is proportional to the ratio of the bursting periods of the uncoupled HCOs, which are no longer equal in the heterogeneous case.

## Iv Coupled Inhibitory CPGs

Restoring the contralateral inhibitory connections, , feed-forwarded from the driving HCO2 to the driven HCO1 enhances the robustness of the bursting patterns of such CPG. Depending on whether it is comprised of homo- or heterogeneous HCOs, the phase-lags between the HCOs can vary, thus giving rise to principally distinct patterns. The unidirectional inhibition gives rise to a polarity in the network. As a result, burst timing of the driven HCO1 has to adjusts itself to that of the driving HCO2 in order for the network to settle into a steady rhythm, if any, with all phases locked. Contralaterality of such inhibition is significant due to the effect it has on timing of interactions when (balancing) ipsilateral excitation is rewired in the CPG schematics as in Fig. 2.

Figure 7 represents the -projections of the phase-lag maps for homogeneous (a) and heterogeneous (c) networks.

In the former case, the majority of the transients tend to a single fixed point at . Still, there is a rudiment of the invariant curve segment, , due to equal convergence rates in the HCOs. These values of the fixed point coordinates are supported by inspection of a few delegated phase-leg progressions plotted against the burst cycle number. Figure 7 yields the frequency count distribution (FCD) of the network states after 100 cycles. The diagram shows a sharp peak at , along with wider peaks (black and blue) at and . The values of the coordinates of the fixed point mean that unidirectional inhibition shifts anti-phase bursting in the driven HCO1 a quarter of the period forward relative to that of the driving HCO2 in the homogeneous network with the same synaptic conductances.

The heterogeneous CPG demonstrates other phase-lags between the HCOs. Recall that in this case the reciprocal inhibitions in HCO1 are halves of those in HCO2: . The corresponding phase-lag map in the -projection given in Fig. 7 shows no indication of an invariant line but the occurrence of a stand-alone stable fixed point. The strong stability of the fixed point is also supported by the fast convergence of transients to the steady states (Fig. 7) after about 10 burst cycles, in contrast to 35 in the homogeneous case. The sharp peaks, at or near for and , and near for in the frequency count distribution of terminal phase-lag is the secondary backup for this assertion. Stated another way, the coordinates, , of this fixed point correspond to the CPG rhythm where both HCOs burst in-phase with the common period.

Next, we will explore how the uni-directional contralateral inhibition and ipsilateral excitation can balance out the CPG dynamics when introduced separately. This should let us identify independent contributions of the synapses of each type to the behavior of the whole CPG.

## V Basic heterogeneous inhibitory and excitatory networks

Intracellular recordings from the four identified interneurons of the swimming CPG of the Melibe have indicated that the phase-lags, (), between the burst initiation in the voltage traces are maintained stably at these values: (). While it is evident that the both HCOs always remain bursting in anti-phase, it is less clear what mechanisms, involving reciprocal inhibition and/or excitation, polarity of wiring etc., are used by 4-neuron networks to preserve the stability of phase-lag between the HCOs. Assuming that the building blocks of such networks remain HCOs formed by anti-phase bursting interneurons, our next step in exploration of such networks, is the examination of functions of contralateral inhibition (circuit in Fig. 8) and ipsilateral excitation (the circuit shown in Fig. 8), and whether they can break or enforce the robustness of activity patterns. In what follows, we explore the dependence of the phase-lag between the HCOs, i.e. , or equivalently, , on the coupling strengths in heterogeneous networks. Panels of Fig. 8 shows how varies as unidirectional feed-forward inhibition from HCO2 onto HCO1, and unidirectional backward excitation from HCO1 onto HCO2 are increased. Since both neural configurations are sub-circuits of the Melibe CPG, let us consider them independently first.

A phase-locked network state turns out to be quantitatively a function of unidirectional inhibition and excitation coupling. When the inhibitory strength, , of the contralateral synapses is increased, the phase-locked state, initially emerging at low small values of , quickly moves to a high value around , see Fig. 8. This means that on average, interneuron 3, first following interneuron 1, becomes delayed by nearly a bursting period as the inhibition between the HCOs substantially increased.

In the case of ipsilateral excitation, increasing works the other way around for synchronization of the interneurons: 1 with 3, and 2 with 4. The evolution of the steady states of the phase-lag, is presented in Fig. 8. It shows that above a threshold , the phase-lag, , shifts down to a steady state, i.e. the driven interneuron 3 (4) follows interneuron 1 (2) after a short delay of -th of the period of the network, and so does HCO2 after HCO1, as a whole.

Given that phase-lags are defined on modulo 1, one can say that in-phase synchrony between the HCOs is due to repulsion in the case of contralateral inhibition, and due to attraction in the excitatory case. A simple calculation (given in appendix) demonstrates that for the network to achieve a robust phase-locked state, the driven HCO has to adjust its period, i.e. either catch up or slow down, in order to match up with that of the driving HCO in unidirectional cases. The effect of increasing synaptic strengths saturates in both cases, after some thresholds are reached. Making the coupling strength five times stronger than the nominal value of the maximal synaptic conductance, has little effect (purple lines in Figs.8 and 8) on the steady state value of . Comparison of two types of coupling (network configurations) suggests that the contralateral inhibition produces a phase-locked state that appears to be the closest to the experimentally observed pattern.

Because there are other, excitatory and electrical connections between the interneurons in the CPG circuitries, in the following sections we will address and identify their roles for predominance and robustness of specific bursting states in 4-cell CPG networks. It is shown in (19) that while the swim CPG of Melibe with strong contralateral inhibitory synapses produces patterns with high values, the swim CPG of a another sea slug Dendronotus, possessing ipsilateral excitatory connections, produces bursting patterns with low values, which agrees with our findings.

### v.1 Modulatory effect of electrical coupling

Electrical coupling, or gap junctions, provide bidirectionally a continuous interaction between interneurons thus affecting synchronization properties of oscillatory neural networks (27). Its magnitude, proportional to the difference between the current values of the membrane potentials, promote in-phase synchronization, in most cases. We introduce electrical coupling (represented by a resistor in the circuitry in Fig. 2) between interneurons 1 and 2 of HCO1 in the form and visa versa, in addition to the contralateral inhibition from HCO2 to HCO1 in the heterogeneous network presented in the previous section. While we refer, in general, to an HCO (HCO1 here) as a pair of interneurons bursting in alternation, a proper gap junctions can overcome the inhibition-caused anti-phase dynamics and synchronize the interneurons to burst as a whole. Prior to that, for values of below a synchronization threshold, the period of HCO1 will gradually vary with an increase of the electrical coupling strength, leading to changes in phase-locked states and transformations of bursting patterns of the CPG, as depicted in Fig. 9.

Figure 9 shows how the phase-lags, , of the CPG change with an increase of the electrical synapse through in HCO1. The gap between interneurons 1 and 2 widens to (toward synchrony at ), as they keep receiving the contralateral inhibition from interneurons 4 and 3, bursting in anti-phase: , as the diagram suggests.

## Vi Range of heterogeneity of the CPG

We have pointed out earlier that weakening reciprocal inhibition between interneurons of one HCO can be equivalent to strengthening reciprocal inhibition in the counterpart. In this section, we examine the range of heterogeneity of the 4-cell network in terms of the misbalance among the synaptic coupling strengths. We will vary the reciprocal inhibitory coupling in HCO1 only, while having those in HCO2 intact along with the contralateral, ipsilateral () and electrical ( ) connections. The following four network configurations, schematically drawn in Fig. 10, are explored in this section. In addition to unidirectional cases, we combine them, in the mixed CPG (Fig. 10), in which the gap junction will next bridge the interneurons of HCO1 (Fig. 10)

A robust phase-locked state of a bursting pattern in a network must persist for a certain range in the high-dimensional parameter space of the coupling weights. Given a large number parameters in a generic 4-cell network with various connections, we need to come up with reduction assumptions to single out one effective control parameter, while other less principle parameters are kept fixed. As such an effective control parameter, we employ the ratio of the inhibitory strengths in individual HCOs. We will start with the case of nearly uncoupled interneurons in HCO1 at small , next the reciprocal inhibition is increased to the nominal value, , and then made 2.5 stronger than . The evolution of the representing phase-lag, , is presented in Fig. 11 for the four network configurations. In the diagram, the purple, black and cyan curves correspond to a representative trajectory converging to an unique attractor as we increase inhibition in HCO1. We note that in all cases in question, both HCOs remain anti-phase bursters, i.e. , in the network, so variations in coupling can only shift the phase-lag between them.

We can conclude from the examination of the convergence tendencies of that a phase-locked state exists in all cases. The widest range of heterogeneity is observed for the feedback configuration in Fig. 11 with the regulatory gap junction. Furthermore, all but one configuration can produce a stable network state with the desired phase-lag .

For the contralaterally inhibitory configuration (I), sketched in Fig. 10, to have , the reciprocal inhibitions within HCO1 must be weaker than in HCO2. Having it equal or stronger, leads to the loss of phase-locked state (transient black and cyan lines in Fig. 11). Fig. 12 presents the outcome of simulations of the network dynamics for longer traces over 350 burst cycles. In this diagram, the only transient (yellow line) for the inhibitory homogeneous CPG (equal coupling) shows no stabilization.

For the configuration (E) with ipsilaterally excitatory synapses between the HCOs (Fig. 10), short-term simulations (45 burst cycles) shows no phase-locking even for the large values of , however continuation simulations over 200 bursting cycles indicate a slow convergence (purple curve in Fig. 12) to the phase-lag .

We can conclude that for the CPG model to maintain robustly and flexibly the desired phase-locking at like that in the Melibe swim CPG, all connections, contralateral inhibitory, ipsilateral excitatory, and electrical, are necessary. These connections provide the feedback loop that widens the range of heterogeneity, within which the CPG network possesses the phase-locked state corresponding to the swim bursting pattern. Another related observation suggests that HCO1 should generate reciprocally inhibition stronger than HCO2 to preserve the locking balance. Finally, addition of relatively strong electrical coupling modulates the phase-locked state (black curve in Fig. 11), that allows the heterogeneous CPG with quite distinct HCOs to maintain the stable phase-lags at .

### vi.1 Reduced maps for the phase-lags between HCOs.

Oscillatory network states and their transformations can be effectively identified and studied through the use of Poincaré return maps. In this section, we use the maps to examine a particular bursting pattern of the inhibitory 4-cell network (in Fig. 10) that corresponds to the in-phase bursting interneurons of the driving HCO2: . In the 2D phase-lag projection, this pattern is associated with the solutions belonging to the main (red) diagonal in Fig. 7. The diagonal is indeed an invariant plane inside the unit cube on which the return for all three phase-lags is defined. In restriction to this plane, the 4-cell network is reduced to the 3-cell one, in which the anti-phase bursting interneurons of HCO1 receive double inhibition during the active phase of the in-phase bursting HCO2. For the reduced network, the return map becomes a two-dimensional map defined on the phase-lags, and .

The dynamics of such a map can be assessed by following forward transients (grey) shown in Fig. 13, whose initial conditions (green dots) are subjected to the synchronization condition . The return map reveals a stable invariant curve wrapping around the unit square (2D torus). In the absence of fixed points, the torus must contain a matching unstable invariant curve too. Because it is unstable and repels forward iterates of the map, we may hypothesize that it wiggles around the unstable in-phase state, , of HCO1. This state is unstable because of breaking perturbations due to periodic forcing originated from HCO2.

In essence, the in-phase bursting HCO2 periodically drives or modulates the phase-lag, , causing the onset of oscillations, or phase jittering, around corresponding to the anti-phase bursting HCO1. This observation lets us define a further reduced 1D map for the discrete evolution of the phase-lag between the reference and any other interneurons of the CPG: , where is the degree of the map. The map for HCO1 is shown in Fig. 13; is chosen because of the slow convergence of transients to an attractor in the weakly coupled case. In this figure, the modulation oscillations of are represented by the closed invariant circle (dark green dots). Such an invariant circle can be made of finite or infinite number of points depending on whether the ratio of the burst periods of the HCOs is rational or not. This circle in the 1D map corresponds to the stable invariant curve wrapping around the 2D torus in Fig. 13. In the case where the interneurons of HCO2 bursts in anti-phase, phase-lag transients converge monotonically to the fixed point at the intersection of the map graph (light green) curve with the -line. The fixed point corresponds to the anti-phase bursting interneurons 1 and 2 of HCO1. An unstable fixed point of the 1D map at the origin (at 1) corresponds to the repelling invariant curve in the 2D map in Fig. 13.

Due to the periodic nature of the patterns in this study, we employ return maps to investigate the ways transients converge to attracting states such as fixed points and invariant circles. Each phase-lag trajectories form a distinct discrete path on the 3D torus (unless it is periodic on an invariant curve). While reducing to 1D return maps is practical in many instances for detecting stable fixed points for phase-locked states of the network, interpretation of solutions of reduced maps can be ill-suited for proper description of high-dimensional dynamics of neural networks. This concerns especially invariant circles and saddle fixed points, which can appear to be stable in restrictions to some invariant subspaces of the 3D maps. As an example, let us discuss the 1D maps shown in Fig. 14. The map in Fig. 14 depicts the transitioning behavior of the forward iterates of (modulo 1) towards a stable fixed point at on the -line. First, the iterates approach from above a phantom at on the -line corresponding to the saddle in the 3D phase space of the full map. Having lingered by, the phase point runs down to re-emerge at the top left corner and next at the right corner of the map. Finally, it converges to the stable fixed point from the right. Figure 14 depicts the behavior of the same sequence in the -projection that tends to the corresponding fixed point at , shifted by a half period. In this projection, the coordinate of the saddle fixed point is .

## Vii Discussion

Many abnormal neurological phenomena are perturbations of normal mechanisms that govern behaviors such as movement. Repetitive motor behaviors are often hypothesized to correspond to rhythmogenesis in small networks of neurons that are able autonomously to generate or continue a variety of activity patterns without further external input. The detailed correspondence across these scales has not yet been made clear in any animal. But there is a growing consensus in the community of neurophysiologists and computational researchers that some basic structural and functional elements are likely shared by CPGs of both invertebrate and vertebrate animals. Before we can study the mechanisms of disorders at the level of individual neurons and CPG circuits in mammals, we therefore first seek to develop better tools and techniques in the context of much simpler animals. However, the ultimate aim of developing our tools and approach to understanding CPGs in lower animals is to make them applicable to studying the governing principles of neurological phenomena in higher animals, and so could potentially assist in treating neurological disorders associated with CPG arrhythmia. Our presentation here is intended as a tutorial guide that demonstrates the effectiveness of our analytical approach that connects exploratory mathematical models to experimental data in the context of known behavioral patterns.

In this pilot study, we focused on a biological CPG that has been linked to a specific swim motion in the sea slug, Melibe. This CPG robustly produces a phase-locked bursting pattern according to recent simultaneous voltage recordings from the four identified neurons. The goal of this modeling study was twofold: first, to identify specific components and their connections in this CPG; second, to get insight into how the connection ‘weights’ in this network lead to robust production of the rhythmic patterns of self-sustained activity. Given the preliminary state of our knowledge about CPGs in the context of whole-animal behavior, we only aim to capture broad rhythmic properties of a small biophysical network without attempting to model each identified Melibe swim CPG neuron in precise detail. Thus, we avoid fitting microscopic details in our model such as a precise constitution of ionic currents, exact shapes and numbers of action potentials per bursts. Instead, we employ a generic Hodgkin-Huxley type model of endogenously bursting neurons that are qualitatively typical for invertebrates. This allows us to concentrate on explore a wide range of network configurations that might be responsible for the specific bursting patterns observed in the Melibe CPG.

We also intentionally construct and present our approach in a fashion that is analogous to the dynamic clamp technique used in neurophysiological experiments. Our technique involves the *dynamic* removal, restoration and variance of (chemical) synaptic connections during simulation, which mimics the experimental techniques of drug-induced synaptic blockade, wash-out, etc. Restoring the chemical synapses during a simulation makes the CPG regain a bursting pattern with temporal characteristics such as phase lags, duty cycles, that depend on the connectivity strengths between the interneurons of the network.

Due to the rhythmic nature of the bursting patterns, we employ Poincaré return maps defined on phases and phase lags between burst initiations in the interneurons. These maps allow us to study quantitative and qualitative properties of the stable rhythms and their corresponding attractor basins. We also exploit conventional knowledge about anti-phase bursting patterns to identify some basic requirements for plausible network configurations. For instance, reciprocal inhibition between a pair of neurons has long been known to produce anti-phase bursts in half center oscillator (HCO) configurations. We rely on a common, standing assumption in current neuroscience that the circuits of motor generation and control are modular in nature. Thus, the present theoretical challenge is how to understand the HCOs as building blocks that must be interconnected to produce single or multiple bursting patterns robustly, and what determines the stability and predominance of these rhythms. Our study is a step towards the dynamical foundation of this theory. We find that our model of the CPG reproduced, quite accurately, the available intracellular recordings from identified interneurons in the Melibe CPG. Furthermore, we find that the Melibe network can be interpreted to consist of two interconnected individual HCOs of two neurons each.

Furthermore, depending on strengths of unidirectional inhibition and excitation, we find that the individual HCOs may have different distributions of phase-locked states. This is a significant observation because, for example, inter-cellular recordings from the identified interneurons of the swim CPG of a similar sea slug, Dendronotus, indicate a phase-locked state that is consistent with our model when it is configured using dominant, ipsilateral excitatory connections from one individual HCO to the other. In addition, as shown in the “Range of heterogeneity” section, the coupling strengths of the reciprocal inhibition within the HCOs have to be balanced in a certain ratio for the whole network to achieve the desired phase-locking.

We dissected the geometric organization of the simulation results into interactions between the building blocks of uncoupled and coupled HCOs. The relations between the phase lags helped us to link network architecture (configuration) with geometric organization of the solutions (model output). For instance, we show that each attractor of the network, whether it is a fixed point or an invariant circle, corresponds to either a phase-locked bursting pattern with distinct phase lags, or else to a bursting pattern with phase lags that vary periodically over the whole network period. As it is unknown, a priori, whether the Melibe swim CPG is multifunctional for given set of parameters, one needs to sweep all possible initial phase distribution to reveal the existence of multiple attractors in the phase space of the corresponding return map. Through the use of decoupled HCOs we were able to explain significant details of the 3D phase portraits such as convergence rates and the occurrence of designated convergence routes to attractors of the phase-lag map.

To study the known robustness of our network, we use an ensemble of computational tools that allow for the reduction of observable voltage dynamics to low-dimensional return maps for phase lags between burst initiations in the interneurons. In particular, we reduce the bursting dynamics of the 4-cell network, represented in full by a 12-dimensional system of ODEs, to a 3-dimensional return map for the phase lags between either the endogenously bursting interneurons or between bursting HCOs. We use a “top-down” approach in which we systematically examine the properties of the phase-lag maps that we abstracted from the 12D system rather than exploring the global dynamics of that full system directly. With this approach, we identified that both contralateral feed-forward inhibition and ipsilateral backward excitation are needed for the network to stabilize the bursting pattern against small perturbations. A certain balance of the synaptic strengths is also required to maintain the phase-locked state within a reasonable range in the parameter space of the CPG network.

Another strong working assumption in this study is that a CPG is composed of (nearly) identical elements — interneurons or bursting HCOs — which are interconnected through chemical synapses and gap junctions of equal conductivity. Due to alterations in the reciprocal wiring, such a homogeneous CPG can be adapted to become dedicated to a single rhythm or multifunctional. We might presume that, through iterative processes of learning and evolution, a real CPG might develop a heterogeneous structure as specific connections become stronger or weaker, so that it can become better adapted to performing specific functions in specific animals. Certainly, we are all aware of examples where, through learning and exercise, mammalian motor systems become “multifunctional” and are able to quickly transition between several dynamic functions on demand: for instance, the diverse swimming styles that have been cultivated by humans, including the in-phase breaststroke and butterfly, and the anti-phase crawl and backstroke. For now, we can only hypothesize that there is a multifunctional, and presumably heterogeneous, CPG network underlying these specific swimming rhythms that determines the phase relationships between rhythmic muscle control signals.

In general, our insights allow us to predict both quantitative and qualitative transformations of the observed patterns whenever the network configurations are altered. The nature of these transformations provides a set of novel hypotheses for biophysical mechanisms about the control and modulation of rhythmic activity. A powerful aspect to our analytical technique is that it does not require knowledge of the equations that model the system. Thus, we believe that have further developed a universal approach to studying both detailed and phenomenological models of bursting networks that is also applicable to a variety of rhythmic biological phenomena beyond motor control.

Even the real Melibe swim CPG is, of course, much more complex than our specific model that is based on the existing, preliminary experimental data. There is great room for improvement in the model by incorporating other biological features into it. There are many open questions that could be addressed by more detailed modeling, such as whether the individual models are natural bursters with distinct duty cycles, or whether they spike tonically and can only become network bursters episodically when under the influence of external drive from other pre-synaptic interneurons in the CPG. This is a challenging question, both phenomenologically and computationally. In future work based upon the framework of this study, we plan to address such questions with more realistic models of inhibitory-excitatory CPG configurations, including ones comprised of three and more HCOs.

## Viii Acknowledgments

This work is funded by NSF grant DMS-1009591, RFFI Grant No. 08-01-00083, MESRF project 14.740.11.0919, as well as GSU Brain and Behaviors Initiative. Research of D. Allen and J. Youker was supported by NSF REU. We would like to thank A. Kelley, J. Schwabedal, R. Clewley, A. Sakurai and P. Katz for helpful discussions and comments. We are also grateful to J. Wojcik and R. Clewley for guidance on the PyDSTool package (37) used for network simulations.

## Ix Appendix

### ix.1 Leech heart interneuron model

This networked model is based on the dynamics of the fast sodium current, , the slow potassium current, , and an ohmic leak current, is given by (7); (15):

(2) |

Here, is the membrane potential, and are the gating variables for sodium channels, which activate (instantaneously) and inactivate, respectively, as the membrane potential depolarizes; is the gating variable for the potassium channel that activates slowly as the membrane potential hyperpolarizes. An applied current, , through the paper unless indicated otherwise. The time constants for the gating variables, maximum conductances and reversal potentials for all the channels and leak current, and the membrane capacitance are set as follows:

The steady state values of the gating variables are given by the following Boltzmann functions:

with V; this bifurcation parameter controls the number of spikes per burst. The currents through fast, non-delayed, chemical synapses are modeled using the fast threshold modulation paradigm as follows (32):

(3) |

here and are the voltages of the post- and the pre-synaptic interneurons; the synaptic threshold V is chosen so that every spike within a burst in the pre-synaptic neuron crosses . This implies that the synaptic current, , is initiated as soon as exceeds the synaptic threshold. The type, inhibitory or excitatory, of the FTM synapse is determined by the level of the reversal potential, , in the post-synaptic neuron. In the inhibitory case, it is set as V so that . In the excitatory case the level of is raised to zero to guarantee that remains, below the reversal potential on average, over the period of the bursting interneuron. We point out specifically that our previous studies of network dynamics revealed that alternative models of fast chemical synapses, using the alpha function and detailed dynamical representation, have no contrast effect on interactions between the coupled interneurons (34).

### ix.2 Phase-lags and HCO periods

The delays between burst initiations in the reference interneuron 1 and other three are given by , and the corresponding times are given by , and recurrent times (network period) by (). The superscripts stand for bursting cycle numbers. Then, the phase-lags are defined as the following:

(4) |

After the phase-locked state is achieved the difference between subsequent phase-lags does not change:

(5) |

The equality, , means that HCO1 maintains a constant period . Then the above conditions can be reduced to those on the following ratios of the periods between both HCOs:

(6) |

Whenever the phase-locked state is achieved, HCO2 has the period of HCO1. On the other hand, if , then the delays and will change proportionally too:

(7) |

### ix.3 Numerical methods

All numerical simulations and phase analysis were performed using the PyDSTool dynamical systems software environment (37). Each sequence of phase lags plotted in Fig. 1 begins from an initial lag , which is the difference in phases measured relative to the recurrence time of cell 1 every time its voltage increases to a threshold mV. marks the beginning of the spiking segment of a burst. As that recurrence time is unknown a priori due to interactions of the cells, we estimate it, up to first order, as a fraction of the period of the synchronous bursting orbit (or that in the individual models) by selecting guess values . The synchronous solution corresponds to . By identifying at the moment when with , we can parameterize this solution by time () or by the phase lag . For weak coupling and small lags, the recurrence time is close to , and . We use the following algorithm to distribute the true initial lags uniformly on a square grid covering the the unit cube (torus), which is the phase space of the phase-lag network.

We initialize the state of cell 1 at from the point of the synchronous solution, and next create the distibution of the initial phase-lagged states in the simulation by suppressing the other cells for durations , respectively. On release, the cells are initialized with the same state . We begin recording the sequence of phase lags between the cells 2–4 and the reference cell 1 on the second cycle after coupling has adjusted the network period away from . In the case of stronger coupling, where the gap between and the first recurrence time for cell 1 widens, we retroactively adjust initial phases using a “shooting” algorithm to make the initial phase lags sufficiently close to uniformly distributed positions on the square grid.

### References

- E. Marder, R. L. Calabrese. Physiol. Rev. 76, 687 (1996);
- J.-M. Goaillard, A.L. Taylor, D. Schulz, D., and E. Marder, Nature Neuroscience, 12, 1424-1430 (2009).
- E. Marder. Introduction to Central Pattern Generators. http://youtu.be/dGRQmndrUzA
- A. Selverston (ed). Model Neural Networks and Behavior, Springer (1985).
- W.B. Kristan, R.L. Calabrese, W.O. Friesen. Prog. Neurobiol. 76, 279 (2005).
- K.L. Briggman, W.B. Kristan. Annu. Rev. Neurosci. 31 (2008).
- A. Shilnikov, R. Gordon, I. Belykh. Chaos 18, 037120 (2008).
- K. Matsuoka. Biol. Cybernetics 56, 345 (1987).
- N. Kopell. Toward a theory of modeling central pattern generators. in Neural Control of Rhythmic Movements in Vertebrates (A.H. Cohen, S. Rossignol, and S. Grillner, eds.) New York, Wiley 369 (1988),
- F. Skinner, N. Kopell, E. Marder, Comput. Neurosci. 1, 69 (1994).
- M. Golubitsky, I. Stewart, P.L. Buono, and J.J. Collins. Nature, 401, 693, (1999).
- I.V. Belykh and A.L. Shilnikov. Phys. Rev. Letters. 101, 078102 (2008).
- W.E. Sherwood, R. Harris-Warrick and J.M. Guckenheimer. J. Comput Neuroscience. 30(2), 323 (2010).
- C. C. Canavier, D. A. Baxter, J. W. Clark, J. H. Byrne. Biol. Cybern. 80, 87 (1999).
- J. Wojcik, R. Clewley and A. Shilnikov. Physics Review E 83, 056209-6 (2011).
- J. Wojcik and A. Shilnikov. http://youtu.be/QNwmp9r279M
- A. Sakurai, J.M. Newcomb, J.L. Lillvis and P.S. Katz. Curr. Biol. 21, 1036 (2011).
- A. Sakurai and P.S. Katz. Soc. Neurosci. Abstr. 37.918.04, (2011).
- J.M. Newcomb, A. Sakurai, J.L. Lillvis, C.A. Gunaratne and P.S. Katz. Proc. Natl. Acad. Sci. USA. 109, 1:10669-76 (2012)
- Y.I. Arshavsky, T.G. Deliagina, G.N. Gamkrelidze, G.N. Orlovsky, Y.V. Panchin, L.B. Popova and O.V. Shupliakov. J Neurophysiol/ 69, 512 (1993).
- Y.I. Arshavsky, T.G. Deliagina, G.N. Orlovsky, Y.V. Panchin, L.B. Popova, R.I. Sadreyev. Ann NY Acad Sci. 16(860), 51 (1998).
- J. Jing, and R. Gillette. J. Neurophysiol. 81(2), 654, (1999)
- W.H. Watson, J.M. Newcomb, S. Thompson. Biol Bull.. 203(2), 152 (2002).
- J.M. Newcomb and W.H. Watson. J. Exp. Biol. 205(3), 397 (2002).
- S. Thompson, W.H. Watson, J. Exp. Biol. 208(7), 1347 (2005).
- E. Marder, Neuron, 76 1 (2012).
- F.K. Skinner, L. Zhang, J.L. Perez Velazquez and P.L. Carlen J. Neurophysiol 81, 1274 (1999).
- A. Shilnikov and G. Cymbalyuk. Phys. Rev. Lett. 94, 048101 (2005).
- P. Chanell, G. Cymbalyuk, and A. Shilnikov. Phys. Rev. Lett. 98, 134101 (2007).
- A. Shilnikov. Nonlinear Dynamics, 68(3), 305 (2012).
- R.J. Calin-Jageman, M.J. Tunstall, B.D. Mensh, P.S. Katz and W.N. Frost. J. Neurophysiol. 98, 2382 (2007).
- N. Kopell and D. Somers. Biol. Cybern. 68, 5 (1993).
- S. Jalil, I. Belykh and A. Shilnikov. Phys. Rev. E. 81, 45201R (2010).
- S. Jalil, I. Belykh and A. Shilnikov. Phys. Rev. E. 85, 36214 (2012).
- C. C. Canavier, R. Butera,R. Dror, J. W. Clark, J. H. Byrne. Biol. Cybern., 68, 1 (1997)
- S. Daun, J.E. Rubin, I.A. Rybak. Comput Neuroscience 27, 3 (2009).
- R. H. Clewley, W. E. Sherwood, M. D. LaMar, and J. M. Guckenheimer. http://pydstool.sourceforge.net (2006).