Asymmetry Factors Shaping Regular and Irregular Bursting Rhythms in Central Pattern Generators

1. Introduction

Central Pattern Generators (CPGs) are neural networks that control muscle function by generating rhythmic patterns (Marder and Bucher, 2001; Grillner, 2003). Many CPGs are built with a non-open network architecture, i.e., a connection topology in which every neuron receives at least one synapse from another member of the CPG (Huerta et al., 2001). The building block in this type of connection architecture is often an oscillator circuit of reciprocally inhibitory neurons (Miller and Selverston, 1982; Selverston, 2010; Sakurai et al., 2014). Another highly relevant property in these neural networks is the presence of intrinsically irregular/chaotic neurons, which typically display rich slow-fast dynamics able to generate bursts of different duration periods, phases and spike temporal structures (Abarbanel et al., 1996; Elson et al., 1998). The reciprocal inhibitory connections between neurons lead to the regularization of the chaotic behavior when the neurons interact with each other in the circuit. The rich intrinsic dynamics provide flexibility and robustness for negotiating rhythms as a function of external inputs, producing the characteristic regular spiking-bursting activity that allows the CPG to control motor movements (Selverston et al., 2000). Most CPGs are typically studied in regular regimes where the activity of individual neurons can be directly related to specific rhythmic motor functions, such as walking, breathing, chewing, etc. However, irregular rhythms can also be observed in control or induced conditions (e.g., see Bartos et al., 1999; Thuma and Hooper, 2003; Nadim et al., 2011; Elices and Varona, 2015; Hooper et al., 2015). Central and sensory feedback transiently alter the CPG pattern which also results in observed irregularity.

Minimal circuits, such as half-center oscillators involving mutually inhibitory neurons are convenient networks to address CPG function. In particular, the concept of a half-center oscillator has been extensively used to study CPG rhythm generation, both in experimental, (e.g., see Miller and Selverston, 1982; Sharp et al., 1996; Yakovenko et al., 2005; Brookings et al., 2012; Sakurai et al., 2014), and theoretical works, (e.g., Nadim et al., 1995; Cymbalyuk et al., 2002; Bem and Rinzel, 2004; Wojcik et al., 2014; Reyes et al., 2015; Doloc-Mihu and Calabrese, 2016). Most of these studies focus in the analysis of alternating regular rhythms, and only a few works address the presence of irregular spiking-bursting activity in mutually inhibitory neurons, (e.g., Varona et al., 2001a; Doloc-Mihu and Calabrese, 2011; Nagornov et al., 2016).

Using conductance-based models and realistic connection topologies inspired by the crustacean pyloric CPG, in this paper we address the study of asymmetry connectivity factors shaping CPG spiking-bursting rhythms. In particular, we show the existence of asymmetric maximal synaptic conductances and time constants that shape regularized and robust alternating spiking-bursting activity in half-center oscillator circuits, and we assess their role in the rhythm configuration. We also discuss the modulation of the regularity by additional gap-junction connections to the half-center oscillator. Finally, we show how a closed-loop protocol can adapt online the synaptic time constants based on the regularity of the burst periods. We argue that this closed-loop interaction is an effective methodology to characterize the coordination properties that arise both from the connection topology and the individual dynamics of the spiking-bursting neurons in these circuits when there is no time to explore the whole parameter space, as in experimental approaches.

2. Materials and Methods

2.1. Neuron Model

For our analysis we use a conductance based model proposed by Komendantov and Kononenko (1996), which describes the conductances of eight membrane currents. The basic membrane potential equation is:

$\begin{array}{cc}\begin{array}{l}-C_{m}dV/dt=I_{Na(TTX)}+I_{K(TEA)}+I_K+I_{Na}+I_{Na}(V)\\ +I_B+I_{Ca}+I_{Ca-Ca}.\end{array}& \left(1\right)\end{array}$

The slow-wave generating mechanism is given by sodium, potassium and chemosensitive currents:

$\begin{array}{cc}{I}_{Na}(V)=g_{Na}^{\ast }(V)(1/(1+exp(-0.2(V+45))))(V-V_{Na});& \left(2\right)\end{array}$

$\begin{array}{cc}{I}_{Na}=g_{Na}^{\ast }(V-V_{Na});& \left(3\right)\end{array}$

$\begin{array}{cc}{I}_K=g_K^{\ast }(V-V_K);& \left(4\right)\end{array}$

$\begin{array}{cc}{I}_B=g_B^{\ast }{m}_B{h}_B(V-V_B);& \left(5\right)\end{array}$

$\begin{array}{cc}d{m}_B/dt=(1/(1+exp(0.4(V+34)))-m_B)/0.05;& \left(6\right)\end{array}$

$\begin{array}{cc}d{h}_B/dt=(1/(1+exp(-0.55(V+43)))-h_B)/1.5;& \left(7\right)\end{array}$

The spike-generating mechanism is described by TTX-sensitive sodium and TEA-sensitive potassium Hodgkin-Huxley type currents:

$\begin{array}{cc}{I}_{Na(TTX)}=g_{Na(TTX)}^{\ast }{m}^{3}h(V-V_{Na});& \left(8\right)\end{array}$

$\begin{array}{cc}{I}_{K(TEA)}=g_{K(TEA)}^{\ast }{n}^4(V-V_K);& \left(9\right)\end{array}$

$\begin{array}{cc}dm/dt=(1/(1+exp(-0.4(V+31)))-m)/0.0005;& \left(10\right)\end{array}$

$\begin{array}{cc}dh/dt=(1/(1+exp(0.25(V+45)))-h)/0.01;& \left(11\right)\end{array}$

$\begin{array}{cc}dn/dt=(1/(1+exp(-0.18(V+25)))-n)/0.015;& \left(12\right)\end{array}$

The calcium transient voltage-dependent current is described by:

$\begin{array}{cc}{I}_{Ca}=g_{Ca}^{\ast }{m}_{Ca}^2(V-V_{Ca});& \left(13\right)\end{array}$

$\begin{array}{cc}d{m}_{Ca}/dt=(1/(1+exp(-0.2V))-m_{Ca})/0.01;& \left(14\right)\end{array}$

The calcium stationary [Ca²⁺]_in inhibited current is given by:

$\begin{array}{cc}\begin{array}{l}{I}_{Ca-Ca}=g_{Ca-Ca}^{\ast }\frac{1}{1+exp(-0.06(V+45))}\\ \frac{1}{1+exp(K_{\beta }([Ca]-\beta ))}(V-V_{Ca});\end{array}& \left(15\right)\end{array}$

$\begin{array}{cc}d[Ca]/dt=\rho -I_{Ca}/2Fv-K_s[Ca];& \left(16\right)\end{array}$

where v = 4πR³/3 is the volume of the cell; [Ca] is ${\left[C{a}^{2+}\right]}_{in}$ (mM), F is Faraday number (F = 96,485 Cmol⁻¹), K_s is the intracellular Ca-uptake rate constant and ρ is the endogenous Ca buffer capacity.

This Hodgkin-Huxley type model displays the dynamical richness observed in the spiking-bursting activity of several neuron types (e.g., see Komarov et al., 2008; Latorre et al., 2013) which underlies some of their rhythm negotiating properties. The parameters used in our simulations are set for the chaotic bursting regime (see Table 1). Figure 1 illustrates the irregular bursting activity of the model in this regime, which resembles the observed behavior in isolated CPG neurons (Abarbanel et al., 1996; Elson et al., 1998).

TABLE 1

Table 1. Parameters for the chaotic bursting regime used in our simulations, (see also appendix in Komendantov and Kononenko, 1996).

FIGURE 1

Figure 1. Chaotic bursting activity of a single Komendantov-Kononenko model neuron. Model parameters are specified in Table 1 (see also Komendantov and Kononenko, 1996). V⁰ = −47.0.

2.2. Network Topologies

Synaptic asymmetries both in strength and duration are known to be present in the mutual inhibition between CPG neurons (Marder and Eisen, 1984). For our study on the asymmetry factors shaping regular and irregular bursting rhythms, we will use three different topologies of mutually inhibitory neurons that we describe below.

2.2.1. Mutually Inhibitory Oscillator Circuit with Symmetric Synaptic Temporal Dynamics

We first consider a half-center oscillator topology, i.e., a minimal network built up with two model neurons connected with mutually inhibitory chemical synapses that have symmetric temporal characteristics but considering the possibility of different maximal conductances in each synapse (see left panel in Figure 2A). The associated synaptic currents have been modeled with a fast graded synapse, a common type of chemical synapse in many CPGs, as follows (Golowasch et al., 1999; Latorre et al., 2002, 2006).

$\begin{array}{cc}{I}_{post}^f=\frac{g_{prepost}^f(V_{post}-E_{Syn})}{1.0+exp(s^f(V^f-V_{pre}))},& \left(17\right)\end{array}$

where $g_{prepost}^f$ is the maximal synaptic conductance of the postsynaptic neuron, one of our control parameters. V_post is the membrane potential of the postsynaptic neuron, E_Syn is the synaptic reversal potential, V^f determines the threshold of the graded synapse, and V_pre is the membrane potential of the presynaptic neuron. The values of the synapse parameters used in our simulations are E_Syn = −65mV, V^f = −49mV, s^f = 0.31mV⁻¹. The mutually inhibitory connections in the circuit lead the neurons to a rhythm negotiation in the form of alternating bursting activity. The activity produced by the circuit can be regular or irregular depending on the value of the maximal conductances of the synapses as it can be observed in the right panels of Figure 2A. In this paper we have considered that the activity is regular if the coefficient of variation of the period C_v during five consecutive bursts is <5%.

FIGURE 2

Figure 2. The three different topologies of mutually inhibitory spiking-bursting neural circuits considered in this study. The activity can be irregular or regular depending on the maximal conductances of the synapses and the rest of parameters that control the asymmetry level in each circuit. (A) mutually inhibitory oscillator circuit with temporal symmetric synapses. Black filled circles represent fast graded chemical synapses, see Equation (17). Parameters: Top panel g₂₁ = 0.005 μS, g₁₂ = 0.044 μS (irregular rhythm); Bottom panel g₂₁ = 0.022 μS, g₁₂ = 0.024 μS (regular rhythm). (B) Mutually inhibitory oscillator circuit with temporal asymmetric synapses. The empty circle represents a slow chemical synapse, see Equations (18–19). Parameters: k₁ = 0.6 ms⁻¹, k₂ = 0.27 ms⁻¹, Top panel g₂₁ = 0.0155 μS, g₁₂ = 0.0255 μS; Bottom panel g₂₁ = 0.0235 μS, g₁₂ = 0.0055 μS. (C) Asymmetric topology with gap-junctions modeled with Eq. (20). Parameters: k₁ = 0.6 ms⁻¹, k₂ = 0.27 ms⁻¹, Top panel g_{21, 31} = 0.05 μS, g_{12, 13} = 0.03 μS; Bottom panel g_{21, 31} = 0.035 μS, g_{12, 13} = 0.055 μS.

2.2.2. Mutually Inhibitory Oscillator Circuit with Asymmetric Synaptic Temporal Dynamics

Temporal asymmetry is present between the LP and the PD neurons in the crustacean pyloric CPG (Marder and Eisen, 1984). For a further characterization of the role of asymmetry shaping CPG spiking-bursting rhythms, one of the fast graded inhibitory synapses in the previous circuit is replaced by a slow graded inhibitory synapse (see right panel in Figure 2B). The slow synaptic current is given in our model by Golowasch et al. (1999) and Latorre et al. (2002, 2006).

$\begin{array}{cc}{I}_{post}^s=g_{prepost}^s{m}_{post}^s(V_{post}-E_{Syn}),& \left(18\right)\end{array}$

where

$\begin{array}{cc}\frac{d{m}_{post}^s}{dt}=\frac{k_1(1.0-m_{post}^s)}{1.0+exp(s^s(V^s-V_{pre}))}-k_2{m}_{post}^s& \left(19\right)\end{array}$

Here $g_{prepost}^s$ is the maximal synaptic conductance of the postsynaptic neuron, k₁ and k₂ are time constants which control the speed and duration of the synaptic current, and V^s determines the threshold of the graded slow synapse. We used V^s = −49 mV and S^s = 1.0 mV⁻¹. Below we will use k₁ and k₂ as additional control parameters to introduce asymmetry in the temporal evolution of the synaptic currents. Examples of regular and irregular alternating bursting activity generated by this type of circuit are shown in the left panels of Figure 2B.

2.2.3. Topology with Gap-Junction Induced Asymmetry

In addition to the asymmetric inhibitory connectivity in the network, we considered a third circuit in which neuron 2 of the previous circuit is replaced by two electrically coupled cells, as it is illustrated in the left panel of Figure 2C. This topology with electrical coupling is also present in several CPGs and in particular in the crustacean stomatogastric ganglion (Selverston and Moulins, 1987). The equation used to model the gap junction was:

$\begin{array}{cc}{I}_{post}^{gj}=g_{prepost}^{gj}(V_{post}-V_{pre}),& \left(20\right)\end{array}$

where the value for the gap-junction conductance used was $g_{prepost}^{gj}=0.0045mS$. Note that this third topology involves two slow inhibitory graded synapses to neuron 1 (see Figure 2C), and one fast graded inhibitory synapse to each of the electrically coupled cells (neurons 2 and 3). Again, the values of the maximal synaptic conductances determine the regular or irregular rhythm in the circuit (see right panels in Figure 2C).

3. Results

3.1. Regularized Activity in Mutually Inhibitory Oscillator Circuits

As a departing point, we have considered a half-center oscillator topology with symmetric fast temporal dynamics in the synapse model (see left panel in Figure 2A). This temporal symmetric circuit generates regular and irregular alternating bursting activity depending on the values of synaptic maximal conductances. Simulations were run to explore the synaptic maximal conductances that resulted in regular or irregular activity according to the criterion of having the coefficient of variation C_v of this activity below a 5% threshold. Figure 3 illustrates the map of conductances g₁₂ and g₂₁ that lead to regularized rhythms (black dots) consisting of alternating bursting between the two neurons. White spaces represent regions were irregular spiking-bursting activity exists. One can observe that this map is nearly symmetric, i.e., regions that lead to regularized spiking-bursting rhythms correspond to values of maximal conductances g₁₂ and g₂₁ that do not differ much beyond 0.02 μS. Thus, for temporal symmetric fast inhibitory synapses, regularization occurs for closely symmetric maximal conductance values.

FIGURE 3

Figure 3. Map of conductances g₁₂ and g₂₁ that lead to regularization of the spiking-bursting activity for the connection topology with fast symmetric temporal dynamics. The incremental step for both g₁₂ and g₂₁ is 0.0025 μS. Black circles represent values that lead to regular activity (C_v < 5%). The symmetric distribution of the circles in the map reflects the balanced temporal dynamics of the fast synapses, which in general allows only moderate differences in the maximal conductances to achieve regularity.

The right panel in Figure 2 shows two examples of regular and irregular alternating bursting activity for representative values of g₁₂ and g₂₁ conductances. Note the symmetric phase relationships of the regular spiking bursting rhythm, which corresponds to the nearly symmetric values of the maximal conductances.

3.2. Regularized Activity in Mutually Inhibitory Oscillator Circuits with Temporal Asymmetry

Here we focus on the analysis of a mutually inhibitory circuit with temporal synaptic asymmetry. In this case, we consider mutual inhibition with a slow synapse in one direction and a fast synapse in the other (see left panel in Figure 2B). This asymmetry is present in half center oscillators that built up CPGs (Marder and Eisen, 1984). Figure 4 represents the map of conductances g₁₂ and g₂₁ that lead to regularized rhythms consisting of alternating spiking-bursting activity between the two neurons. It is important to emphasize that in this case asymmetry in the maximal conductances coexists with the temporal asymmetry of the slow and fast synapses in the circuit. Note the reduced size and the sparser distribution of maximal conductances that lead to regularized spiking-bursting rhythms, which is caused by the synaptic imbalance introduced by the temporal asymmetry. Figure 2B illustrates two examples of regular and irregular spiking-bursting activity, respectively, for representative values of g₁₂ and g₂₁ conductances. One can see the asymmetric phases of the regular rhythm as compared to the one depicted in Figure 2A.

FIGURE 4

Figure 4. Map of conductances g₁₂ and g₂₁ of the temporal asymmetric mutually inhibitory circuit that lead to regularization of the spiking-bursting activity. The incremental step for both g₁₂ and g₂₁ is 0.0025 μS. Black circles represent values that lead to regular activity (C_v < 5%). The map shows a reduced number of configurations of parameters that leads to regular activity as a consequence of the temporal asymmetry in the connection between the neurons. Parameters: k₁ = 0.6 ms⁻¹, k₂ = 0.27 ms⁻¹. Red squares correspond to specific values of g₁₂ and g₂₁ that will be used for subsequent analysis of k₁ vs k₂ maps.

Next, we chose representative values of the maximal conductances in the temporal asymmetric circuit (indicated by red squares in Figure 4) and, using the same methodology, we explored the parameter space of the time constants k₁ and k₂, which control the temporal aspects of the slow synapse. Figure 5 shows this analysis for fixed g₂₁ = 0.03μS and g₁₂ = 0.0125 μS, corresponding to a regularized activity regime in Figure 4. Note the presence of k₂ bands where regularity occurs for large regions of k₁ values for this particular selection of maximal conductances. The size of these bands, and thus the size of the regions with regular activity, depend on the values of the maximal conductances of the mutually inhibitory circuit.

FIGURE 5

Figure 5. Map of slow synaptic time constants k₁ and k₂ of the temporal asymmetric mutually inhibitory oscillator circuit which lead to regularization of the spiking-bursting activity. Parameters: g₂₁ = 0.03 μS, g₁₂ = 0.0125 μS (see Figure 4 red squares). The incremental step for k₁ and k₂ is 0.025 ms⁻¹ and 0.005 ms⁻¹, respectively. Black circles represent values that lead to regular activity (C_v < 5%).

Figure 6 shows the k₁ vs k₂ map for fixed g₂₁ = 0.065 and g₁₂ = 0.005 μS. Note that such conductance values also correspond to regular spiking-bursting activity in the map depicted in Figure 4. In spite of the large maximal conductance difference, the resulting combination of conductance and temporal asymmetries in the mutual inhibition lead to a broad region of regular spiking-bursting activity. Irregular regimes are mostly in the region defined by k₂ < 0.3

FIGURE 6

Figure 6. Map of slow synaptic time constants k₁ and k₂ of the temporal asymmetric mutually inhibitory oscillator circuit which lead to regularization of the spiking-bursting activity. Parameters: g₂₁ = 0.065 μS, g₁₂ = 0.005 μS (see Figure 4 red squares). The incremental step for k₁ and k₂ is 0.025 ms⁻¹ and 0.005 ms⁻¹, respectively. Black circles represent values that lead to regular activity (C_v < 5%).

When we select from Figure 4 a set of maximal conductances that correspond to irregular spiking-bursting activity (e.g., g₂₁ = 0.065 and g₁₂ = 0.015 μS), we obtain a k₁ vs k₂ map with reduced regularity bands (see Figure 7). Nevertheless, one can observe that the temporal asymmetry can compensate the conductance unbalance to achieve regular regimes.

FIGURE 7

Figure 7. Map of slow synaptic time constants k₁ and k₂ of the temporal asymmetric mutually inhibitory oscillator circuit which lead to regularization of the spiking-bursting activity. Parameters: g₂₁ = 0.065 μS, g₁₂ = 0.015 μS (see Figure 4 red squares). The incremental step for k₁ and k₂ is 0.025 ms⁻¹ and 0.005 ms⁻¹, respectively. Black circles represent values that lead to regular activity (C_v < 5%).

3.3. Regularized Activity in Mutually Inhibitory Oscillator Circuits with Gap-Junction Induced Asymmetry

Finally, in addition to the asymmetric connectivity in the network, we have considered another source of asymmetry common in many CPG's (Selverston and Moulins, 1987) by replacing one of the neurons in the former circuit with two electrically coupled cells (see left panel in Figure 2C). The gap-junction synchronizes the activity of neurons 2 and 3. It is important to note that in this configuration there are two slow and two fast inhibitory graded synapses. Note also the phase asymmetry of the regular regimes in this configuration. To analyze the contribution of this topological asymmetry to the generation of regular and irregular rhythms, we have fixed the synaptic time constants to k₁ = 0.6 ms⁻¹, k₂ = 0.27 ms⁻¹, and the gap-junction conductance $g_{23,32}^{gj}$ = 0.0045 μS.

Figure 8 shows the map of conductances g₁₂ and g₂₁ that lead to regularized alternating bursting activity in this configuration. One can observe that the regularity regions are larger than in the corresponding case for the circuit with only temporal synaptic asymmetry (cf. Figure 4). This result hints that the gap-junction connectivity, with its associated dynamical inertia, could contribute to enlarge the regions of regular spiking-bursting activity in living half-center oscillators.

FIGURE 8

Figure 8. Map of conductances g₁₂ and g₂₁ of the circuits with gap-junction induced asymmetry that lead to regularization of the spiking-bursting activity where the incremental step for both g₁₂ and g₂₁ is 0.0025 μS. Black circles represent values that lead to regular activity (C_v < 5%). Parameters: k₁ = 0.6 ms⁻¹, k₂ = 0.27 ms⁻¹, $g_{23,32}^{gj}$ = 0.0045 μ S.

4. Closed-loop Exploration for Regular Bursting Rhythms

In the previous sections we have presented a modeling study in minimal mutually inhibitory circuits of the asymmetry factors that contribute to shape regular and irregular spiking-bursting rhythms. How could the hypotheses drawn in this paper be tested experimentally? Dynamic clamp can be used to modify maximal conductances and synaptic temporal characteristics by building hybrid circuits composed of living neurons and artificial synapses (Sharp et al., 1993; Destexhe and Bal, 2009). As the time to perform such experiments is restricted and the number of configurations that can be tested is typically small, we propose here a close-loop exploration method to allow such study. The closed-loop algorithm can control the dynamics by changing a set of parameters of the network to achieve a specific goal (e.g., the regularization of the bursting activity). These set of parameters are updated in every iteration according to a designed rule (which should be simple due to time restrictions) until the goal is achieved.

We validate the method in the context of the half-center oscillator model with temporal asymmetric synapses and g₂₁ = 0.065 μS, g₁₂ = 0.005 μS. In this case, the closed-loop algorithm updates the value of the synaptic time constants k_{1, 2} according to the difference between the current period P₁ and the average of the periods of the last five bursts 〈P₁〉 trying to find a suitable new set of values that lead to regular activity C_v < 5%, which is the goal of this closed-loop exploration. The departing points are different sets of configurations of the synaptic time constants $k_1^0$ and $k_2^0$ that drive the system into alternating irregular spiking-bursting activity. Every time a new burst n is generated in neuron 1, the algorithm checks if the coefficient of variation of the period calculated for the last 5 bursts C_v is below the established regularity threshold 5%. Then, the value of the synaptic time constants are updated as follows:

$k_{1,2}^n=k_{1,2}^{n-1}-\alpha ·\left(\langle P_1\rangle -P_1^n\right)$

where α is the rate for the synaptic time constant change, 〈P₁〉 is the average period of neuron 1 calculated with the last five bursts and $P_1^n$ is the current period. In this study we set α = 0.01. This value and the monitoring of the last five bursts provide a good balance between a reasonable online measurement of the regularity and a fast convergence to the closed-loop goal.

Figure 9 shows the evolution of three sets of departing time constants $k_1^0$ and $k_2^0$ during the closed-loop protocol displayed over a map of slow synaptic time constants that lead to regularized rhythms (gray circles). In these three examples the departing values are located in empty regions in the maps (C_v > 5%) and at the end of the protocol the new set of values $k_1^f$ and $k_2^f$ are located within the gray circles area (C_v < 5%). Noting the color code, one can observe that the changes in $k_{1,2}^n$ can become quite large during the first few interactions since this change is proportional to the difference between the current period and the average of the periods of the last five bursts. The protocol is able to find the regular spiking-bursting regimes in only a few iterations. The departing irregularity and the achieved regularized activity of the circuit after the closed-loop protocol are shown in Figure 10.

FIGURE 9

Figure 9. Evolution of the synaptic time constants k₁ and k₂ of the low synapse in a half-center oscillator during closed-loop exploration for regular bursting rhythms. The figure shows the evolution of the values of k₁ and k₂ in time (color bar in s) for three different departing values (coded in black): squares points, $k_1^0$ = 1.00 ms⁻¹, $k_2^0$ = 0.63 ms⁻¹; triangle points, $k_1^0$ = 0.70 ms⁻¹, $k_2^0$ = 0.25 ms⁻¹; circle points, $k_1^0$ = 1.80 ms⁻¹, $k_2^0$ = 0.27 ms⁻¹. Final values after the closed-loop exploration (coded in yellow): squares points, $k_1^f$ = 0.92 ms⁻¹, $k_2^f$ = 0.55 ms⁻¹; triangle points, $k_1^f$ = 0.82 ms⁻¹, $k_2^f$ = 0.37 ms⁻¹; circle points, $k_1^f$ = 1.89 ms⁻¹, $k_2^f$ = 0.35 ms⁻¹. Gray circles in the background represent the map of slow synaptic time constants k₁ and k₂ that lead to regularization of the spiking-bursting activity. This map was generated as explained in the previous sections and is depicted for the validation of the closed-loop search. The incremental step for k₁ and k₂ is 0.05 ms⁻¹ and 0.01 ms⁻¹, respectively. The protocol searchers for a suitable value of both synaptic time constants which drives the system into a regular bursting activity, C_v < 5% (gray circle points). Parameters: g₂₁ = 0.065 μS, g₁₂ = 0.005 μS.

FIGURE 10

Figure 10. Departing irregular regimes and closed-loop exploration for regular activity for the three cases discussed in Figure 9 in a circuit with synaptic temporal asymmetry. Left panels show the departing irregular activity for the following slow synaptic time constants: $k_1^0$ = 1.00 ms⁻¹, $k_2^0$ = 0.63 ms⁻¹, $k_1^0$ = 0.70 ms⁻¹, $k_2^0$ = 0.25 ms⁻¹ and $k_1^0$ = 1.80 ms⁻¹, $k_2^0$ = 0.27 ms⁻¹. Right panels show the regularization obtained through the closed-loop (cf. evolution of time constants in Figure 9). Note that the protocol is able to find the regular regimes in a very short time. Maximal conductance parameters: g₂₁ = 0.065 μS, g₁₂ = 0.005 μS.

This close-loop protocol allows different configurations and optimizations. The experimenter can choose to modify time constants k₁ and k₂ simultaneously or one at a time, for example only k₂ to unveil the regularity bands shown in Figures 5, 7. The protocol can also include a nonlinear dependence on the difference between the current and the averaged bursting period, as well as other search methods (e.g. gradient descent, stochastic search) and alternative performance measurements to evaluate online the regularization goal.

5. Discussion

Using realistic conductance-based models of half-center oscillator based CPG circuits, in this paper we have addressed the study of synaptic asymmetry to shape robust alternating spiking-bursting rhythms. In particular, we focused on the analysis of asymmetry of synaptic maximal conductances and time constants in the reciprocal inhibition of half center oscillators. We have mapped the regimes of regular and irregular coordinated rhythms as a function of these parameters. We have shown that asymmetry both in the maximal conductances and in the temporal dynamics of mutually inhibitory neurons, including modulation by gap-junction connectivity, can synergistically contribute to shape large regions of regular spiking-bursting regimes in central pattern generator circuits. Regular rhythms resulting from realistic asymmetric synaptic configurations display specific phase relationships that reflect the balance among the distinct sources of asymmetry.

Irregularity regimes are typically disregarded both in experimental and theoretical CPG research. One reason for this is that regular rhythms, as recorded in experimental setups, are more easily associated with observable rhythmic CPG motor functions. Several theoretical studies have shown that regular regimes could be more efficient for specific motor tasks (Huerta et al., 2000, 2001; Stiesberg et al., 2007). However, coordinated irregularity is also present in living CPGs under normal and pathological circumstances and might mediate key aspects of the rhythm negotiation in these circuits. Recent brute force approaches to map the parameter space in half center oscillators models inspired by the leech heartbeat CPG have also pointed out the presence of irregular regimes (Doloc-Mihu and Calabrese, 2011, 2016). In this paper, we have shown that the asymmetric synaptic parameter space for the existence of regular and irregular coordinated rhythms in a conductance-based half-center oscillator is large. Previous modeling efforts that did not explore synaptic asymmetry hinted that intrinsic neuron irregularity easily disappeared under mutual inhibition in CPG half-center oscillator circuits, (e.g., see Varona et al., 2001a,b). Here we have seen how half-center oscillations can be modulated and coordinated by asymmetrical factors in the mutual inhibition of its constituent neurons. The study of irregular regimes in the case of realistic asymmetric synapses can shed further light to understand the balance between robustness and flexibility in central pattern generator circuits.

There are known asymmetries in the connectivity of CPG half-center oscillators. However, their role in shaping the circuit rhythm has not been explored in detail. In this paper, we have shown that a closed-loop protocol adapting online the synaptic parameters under a regularization goal can be an effective methodology to map and characterize the coordination properties that arise both from the connection topology, including their asymmetry, and the individual neuronal dynamics in these circuits. We believe that customized versions of this protocol can be used in experimental setups, where there is no time to explore the whole parameter space, to address the hypotheses drawn by the discussed model.

Author Contributions

IE and PV designed the study. IE performed the simulations. IE and PV analyzed the data and wrote the paper. All authors approved the final version for publication.

Funding

We acknowledge support from MINECO/FEDER DPI2015-65833-P and TIN-2012-30883 (http://www.mineco.gob.es/) and ONRG grant N62909-14-1-N279.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Abarbanel, H. D. I, Huerta, R., Rabinovich, M. I, Rulkov, N. F., Rowat, P. F., and Selverston, A. I (1996). Synchronized action of synaptically coupled chaotic model neurons. Neural Comput. 8, 1567–1602. doi: 10.1162/neco.1996.8.8.1567

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartos, M., Manor, Y., Nadim, F., Marder, E., and Nusbaum, M. P. (1999). Coordination of fast and slow rhythmic neuronal circuits. J. Neurosci. 19, 6650–6660.

PubMed Abstract | Google Scholar

Bem, T., and Rinzel, J. (2004). Short duty cycle destabilizes a half-center oscillator, but gap junctions can restabilize the anti-phase pattern. J. Neurophysiol. 91, 693–703. doi: 10.1152/jn.00783.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookings, T., Grashow, R., and Marder, E. (2012). Statistics of Neuronal Identification with Open- and Closed-Loop measures of intrinsic excitability. Front. Neural Circ. 6:19. doi: 10.3389/fncir.2012.00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Cymbalyuk, G. S., Gaudry, Q., Masino, M. A., and Calabrese, R. L. (2002). Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J. Neurosci. 22, 10580–10592. doi: 10.3410/f.1011296.179117

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., and Bal, T. (eds.). (2009). Dynamic-Clamp: From Principles to Applications. New York, NY: Springer.

Google Scholar

Doloc-Mihu, A., and Calabrese, R. L. (2011). A database of computational models of a half-center oscillator for analyzing how neuronal parameters influence network activity. J. Biol. Phys. 37, 263–283. doi: 10.1007/s10867-011-9215-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Doloc-Mihu, A., and Calabrese, R. L. (2016). Analysis of family structures reveals robustness or sensitivity of bursting activity to parameter variations in a Half-Center Oscillator (HCO) Model. eNeuro 3:e0015–16. doi: 10.1523/ENEURO.0015-16.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Elices, I., and Varona, P. (2015). Closed-loop control of a minimal central pattern generator network. Neurocomputing 170, 55–62. doi: 10.1016/j.neucom.2015.04.097

CrossRef Full Text | Google Scholar

Elson, R. C., Selverston, A. I, Huerta, R., Rulkov, N. F., Rabinovich, M. I, and Abarbanel, H. D. I (1998). Synchronous behavior of two coupled biological neurons. Phys. Rev. Lett. 81, 5692–5695. doi: 10.1103/PhysRevLett.81.5692

CrossRef Full Text | Google Scholar

Golowasch, J., Casey, M., Abbott, L. F., and Marder, E. (1999). Network stability from activity-dependent regulation of neuronal conductances. Neural Comput. 11, 1079–1096. doi: 10.1162/089976699300016359

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S. (2003). The motor infrastructure: from ion channels to neuronal networks. Nat. Rev. Neurosci. 4, 573–586. doi: 10.1038/nrn1137

PubMed Abstract | CrossRef Full Text | Google Scholar

Hooper, R. M., Tikidji-Hamburyan, R. A., Canavier, C. C., and Prinz, A. A. (2015). Feedback control of variability in the cycle period of a central pattern generator. J. Neurophysiol. 114, 2741–2752. doi: 10.1152/jn.00365.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Huerta, R., Sánchez-Montañès, M. A., Corbacho, F., and Sigüenza, J. A. (2000). A Central pattern generator to control a pyloric-based system. Biol. Cybern. 82, 85–94. doi: 10.1007/PL00007963

PubMed Abstract | CrossRef Full Text | Google Scholar

Huerta, R., Varona, P., Rabinovich, M. I, and Abarbanel, H. D. I (2001). Topology selection by chaotic neurons of a pyloric central pattern generator. Biol. Cybern. 84, L1–L8. doi: 10.1007/pl00007976

PubMed Abstract | CrossRef Full Text | Google Scholar

Komarov, M. A., Osipov, G. V, and Suykens, J. A. K. (2008). Variety of synchronous regimes in neuronal ensembles. Chaos 18:037121. doi: 10.1063/1.2959340

PubMed Abstract | CrossRef Full Text | Google Scholar

Komendantov, A. O., and Kononenko, N. I (1996). Deterministic chaos in mathematical model of pacemaker activity in bursting neurons of snail, helix pomatia. J. Theor. Biol. 183, 219–230. doi: 10.1006/jtbi.1996.0215

PubMed Abstract | CrossRef Full Text | Google Scholar

Latorre, R., Levi, R., and Varona, P. (2013). Transformation of context-dependent sensory dynamics into motor behavior. PLoS Computat. Biol. 9:e1002908. doi: 10.1371/journal.pcbi.1002908

PubMed Abstract | CrossRef Full Text | Google Scholar

Latorre, R., Rodríguez, F. B., and Varona, P. (2002). Characterization of triphasic rhythms in central pattern generators (I): interspike interval analysis. Lect. Notes Comput. Sci. 2415, 167–173. doi: 10.1007/3-540-46084-5-28

CrossRef Full Text | Google Scholar

Latorre, R., Rodríguez, F. B., and Varona, P. (2006). Neural signatures: multiple coding in spiking-bursting cells. Biol Cybern 95, 169–183. doi: 10.1007/s00422-006-0077-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Marder, E. Bucher, D. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, R986–R996. doi: 10.1016/S0960-9822(01)00581-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Marder, E. Eisen, J. S. (1984). Transmitter identification of pyloric neurons: electrically coupled neurons use different transmitters. J. Neurophys. 51, 1345–1361.

PubMed Abstract | Google Scholar

Miller, J. P., and Selverston, A. I (1982). Mechanisms underlying pattern generation in lobster stomatogastric ganglion as determined by selective inactivation of identified neurons. IV. Network properties of pyloric system. J. Neurophysiol. 48, 1416–1432.

PubMed Abstract | Google Scholar

Nadim, F., Olsen, O. H., De Schutter, E., and Calabrese, R. L. (1995). Modeling the leech heartbeat elemental oscillator. I. Interactions of intrinsic and synaptic currents. J. Computat. Neurosci. 2, 215–235. doi: 10.1007/BF00961435

PubMed Abstract | CrossRef Full Text | Google Scholar

Nadim, F., Zhao, S., Zhou, L., and Bose, A. (2011). Inhibitory feedback promotes stability in an oscillatory network. J. Neural Eng. 8:065001. doi: 10.1088/1741-2560/8/6/065001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagornov, R., Osipov, G., Komarov, M., Pikovsky, A., and Shilnikov, A. (2016). Mixed-mode synchronization between two inhibitory neurons with post-inhibitory rebound. Commun. Nonlin. Sci. Numerical Simul. 36, 175–191. doi: 10.1016/j.cnsns.2015.11.024

CrossRef Full Text | Google Scholar

Reyes, M. B., Carelli, P. V, Sartorelli, J. C., and Pinto, R. D. (2015). A modeling approach on why simple central pattern generators are built of irregular neurons. PLoS ONE 10:e0120314. doi: 10.1371/journal.pone.0120314

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakurai, A., Gunaratne, C. a., and Katz, P. S. (2014). Two interconnected kernels of reciprocally inhibitory interneurons underlie alternating left-right swim motor pattern generation in the mollusc Melibe leonina. J. Neurophysiol. 112, 1317–1328. doi: 10.1152/jn.00261.2014

CrossRef Full Text | Google Scholar

Selverston, A. I (2010). Invertebrate central pattern generator circuits. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 2329–2345. doi: 10.1098/rstb.2009.0270

PubMed Abstract | CrossRef Full Text | Google Scholar

Selverston, A. I, and Moulins, M. (eds.). (1987). The Crustacean Stomatogastric System: a Model for the Study of Central Nervous System. Berlin; Heidelberg; New York, NY; London; Paris; Tokyo: Springer-Verlag.

Google Scholar

Selverston, A. I, Rabinovich, M. I, Abarbanel, H. D. I, Elson, R., Szücs, A., Pinto, R. D., et al. (2000). Reliable circuits from irregular neurons: a dynamical approach to understanding central pattern generators. J. Physiol. Paris 94, 357–374. doi: 10.1016/S0928-4257(00)01101-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharp, A. A., O'Neil, M. B., Abbott, L. F., and Marder, E. (1993). The dynamic clamp: artificial conductances in biological neurons. Trends Neurosci. 16:389. doi: 10.1016/0166-2236(93)90004-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharp, A. A., Skinner, F. K., and Marder, E. (1996). Mechanisms of oscillation in dynamic clamp constructed two-cell half-center circuits. J. Neurophysiol. 76, 867–883.

PubMed Abstract | Google Scholar

Stiesberg, G. R., Reyes, M. B., Varona, P., Pinto, R. D., and Huerta, R. (2007). Connection topology selection in central pattern generators by maximizing the gain of information. Neural Comput. 19, 974–993. doi: 10.1162/neco.2007.19.4.974

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuma, J. B., and Hooper, S. L. (2003). Quantification of cardiac sac network effects on a movement-related parameter of pyloric network output in the lobster. J. Neurophysiol. 89, 745–753. doi: 10.1152/jn.00631.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Varona, P., Torres, J. J., Abarbanel, H. D. I, Rabinovich, M. I, and Elson, R. C. (2001a). Dynamics of two electrically coupled chaotic neurons: experimental observations and model analysis. Biol. Cybern. 84, 91–101. doi: 10.1007/s004220000198

PubMed Abstract | CrossRef Full Text | Google Scholar

Varona, P., Torres, J. J., Huerta, R., Abarbanel, H. D. I, and Rabinovich, M. I (2001b). Regularization mechanisms of spiking-bursting neurons. Neural Networks 14, 865–875. doi: 10.1016/S0893-6080(01)00046-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wojcik, J., Schwabedal, J., Clewley, R., and Shilnikov, A. L. (2014). Key bifurcations of bursting polyrhythms in 3-cell central pattern generators. PLoS ONE 9:e92918. doi: 10.1371/journal.pone.0092918

PubMed Abstract | CrossRef Full Text | Google Scholar

Yakovenko, S., McCrea, D. A., Stecina, K., and Prochazka, A. (2005). Control of locomotor cycle durations. J. Neurophysiol. 94, 1057–1065. doi: 10.1152/jn.00991.2004

PubMed Abstract | CrossRef Full Text | Google Scholar