HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 95/3/1556    most recent 01157.2005v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Boothe, D. L. Articles by Troyer, T. W. Search for Related Content PubMed PubMed Citation Articles by Boothe, D. L. Articles by Troyer, T. W.

Temporal Correlations in Stochastic Models of Double Bursting During Simulated Locomotion

¹Neurosciences and Cognitive Sciences Program, ²Department of Biology, ³Institute for Systems Research, and ⁴Department of Psychology, University of Maryland, College Park, Maryland

Submitted 2 November 2005; accepted in final form 8 December 2005

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The output of the spinal central pattern generator for locomotion falls into two broad categories: alternation between antagonistic muscles and double bursting within muscles acting on multiple joints. We first model an alternating half-center and then present two different models of double bursting. The first double-bursting model consists of a central clock with an explicit one-to-one mapping between interneuron activity and model output. The second double-bursting model consists of a half-center with an added feedback neuron. Models are built using rate-coded leaky integrator neurons with slow self-inhibition. Structure-function relationships are explored by the addition of noise. The interaction of noise with the dynamics of each network creates a unique pattern of correlation between phases of the simulated cycle. The effects of noise can be explained by perturbation of deterministic versions of the networks. Three basic results were obtained: slow self-inhibitory currents lead to correlations between parts of the step cycle that are separated in time and network relative; model outputs are most sensitive to perturbations presented just before competitive switches in network activity, and clock-like models possess substantial symmetries within the correlation structure of burst durations, whereas the correlation structure of feedback models are asymmetric. Our models suggest that variability in burst length durations can be analyzed to make inferences about the structure of the spinal networks for locomotion. In particular, correlation patterns within double-bursting outputs may yield important clues regarding the interaction between more central, clock-like networks and feedback from more peripheral interneurons.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The dominant behavior of electroneurograms (ENGs) recorded in the mammalian spinal cord during locomotion is an alternating output that occurs during walking. Generally, flexors are active during the swing phase of locomotion when the animal lifts its limb off the ground, and extensors are active during stance when the animal makes contact with the ground and pushes itself forward. Many muscles display more complicated behaviors such as double bursting, in which a flexor is active during both the swing and the stance phases of locomotion (Grillner 1979). Within intact animals, the period of co-activation between flexor activity and extensor activity is thought to stabilize the limb on contact with the ground and can be used to create additional force during extension (Grillner 1975). The addition of a second flexion burst to the familiar alternation between flexion and extension results in an output pattern with four parts (Forssberg 1979): flexor on alone (F); extensor on alone during the beginning part of extension (E1); a second flexion phase cotemporaneous with extensor activity (E2); and extensor on alone during the later part of extension (E3).

Double bursting raises a key theoretical issue about the structure of the spinal central pattern generator (CPG) for locomotion: are the two bursts of flexion controlled by independent or overlapping sets of neurons within the CPG? According to one hypothesis, the CPG is a set of neurons working much like a clock whose timing is set centrally (Burke 1999; Burke et al. 2001). The basic output pattern of this clock is then fine tuned by the many interneurons of the spinal cord, but the basic pattern of activity within the clock itself remains unchanged. Alternatively, the interaction between more centrally located interneurons and feedback from interneurons closer to the periphery may create the proper output (Edgley and Jankowska 1987; Shefchyk et al. 1990). Under this hypothesis, some spinal premotor interneurons are part of the CPG and contribute to pattern formation even in the absence of sensory input.

Here we investigate three network models of the spinal CPG for locomotion built using mutually inhibitory leaky integrator neurons with slow self-inhibition. We first study the dynamic properties of these networks using a simple half-center model with two mutually inhibitory interneurons. We then extend the half-center to consider networks that embody the two hypotheses above for producing double-bursting outputs. The closest point of experimental contact with our modeling work is the data obtained from fictive walking in which the spinal cord is isolated by removing descending motor commands and inputs from the periphery. Although fictive walking is often quite variable and sometimes nonfunctional, there are bouts of fictive locomotion where the CPG output is very similar to the output of the intact system (Grillner 1975; Grillner and Zangger 1984).

To investigate structure-function relationships in our models, we injected white-noise currents into the model interneurons and studied the resulting patterns of variability in motor neuron outputs. In particular, we measured patterns of variance and covariance in the durations of different parts of the step cycle. To further analyze the models, we performed a series of microstimulation experiments that allows the effect of a single perturbation to be traced through time. Results from theseexperiments can then be reintegrated to understand the results from the noise-driven simulations.

Overall, we find three basic results. First, we find that the slow dynamics of self-inhibition can cause correlations that occur over more than a single cycle. For example, in one network model, the duration of the flexion burst F is strongly correlated with the last phase of extension E3 but is nearly independent of the intervening phases E1 and E2. Second, the microstimulation studies show that the networks are much more sensitive to perturbation at some time points than others. For the current class of models, this sensitivity is generally greatest during the period just before transitions between parts of the step cycle. Finally, we find that the underlying symmetry in the clock-based model is reflected in the variations of the durations of different phases of the cycle. In particular, all phases of the step cycle have similar variance, and the correlation between two phases is determined chiefly by the number of intervening phases of the step cycle rather than the identity of either phase. This pattern is not seen within the model based on feedback from interneurons closer to the periphery where correlations between bursts are specific to interactions between the underlying inhibitory interneurons. Our modeling results suggest that measurements of the noise-driven correlation structure may provide insight into the structure and function of the underlying spinal CPG for locomotion.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The main aim of the models presented here was to explore basic theoretical issues related to the production of double-bursting outputs during locomotion. Given our limited knowledge concerning the neuronal properties and connectivity of mammalian pattern generating circuits, we have chosen to explore basic structure-function relationships in simplified models of the spinal CPG (for current debates see Eide et al. 1999; Huang et al. 2000). The connectivity and structure of the networks are shown in Fig. 1. The connectivity of the half-center model is similar to that proposed for the lamprey spinal CPG (Grillner 1975). The connectivity of the clock model has been studied previously within the context of gait creation and symmetry (Golubitsky et al. 1999), and the connectivity of the feedback model is similar to the function proposed for some sensory interneurons of the spinal cord (Edgely and Jankowska 1987).

View larger version (21K): [in this window] [in a new window]   FIG. 1. Model architectures. Excitatory/inhibitory connections are marked with a line/filled circle. Excitatory interneurons form the model output: Ef, to flexors; Ee to extensors. A: half-center model. Mutual inhibition biased by slow self-inhibition (not shown) causes alternation of the 2 central interneurons (Ie, If) and hence alternating flexion and extension. B: clock model. Mutual inhibition with asymmetric weights causes serial activation of the 4 central interneurons, I1 to I2 to I3 to I4. When I3 is active, neither Ef nor Ee is inhibited, creating double-bursting output. C: feedback model consisting of half-center plus an additional "override" neuron Io. Onset of extension activates Io, inhibiting If and Ie simultaneously, creating double-bursting output.

We used the well-studied lamprey system as a starting point for our networks (Jung et al. 1996; Williams 1992). Lamprey neurons are known to exhibit bursting behavior in which spikes ride on a wave of depolarization. This wave of depolarization is well correlated with the output of the motor neurons as recorded by ENGs (Grillner 1975, 1979, 1999; Kiehn et al. 1996). Therefore we used simple rate-coded neurons where the activity in a single-model neuron represents the summed bursting activity of a population of functionally related neurons (Boothe and Cohen 2001, 2002; Jung et al. 1996). Model neuron i is governed by two internal variables: a normalized membrane voltage variable V_i that takes values between –1 and 1, and a slow self-inhibitory conductance g_diD_i, where D_i ranges between 0 and 1. Voltage dynamics are based on the following differential equation

(1)

The terms on the right hand side of Eq. 1 represent the following currents: a leak current with conductance g_ri and reversal potential V = 0; a tonically active excitatory current with conductance g_ti and reversal potential V = 1; synaptic input from other neurons as represented by the sigmoid firing rate function h(V_j) times a maximal synaptic conductance g_syn^ij and reversal potential v_syn (v_syn = 1 for excitation, v_syn = –1 for inhibition); the slow self-inhibitory current with conductance g_diD_i and reversal potential V = –1; and a perturbation current I_i. Perturbation currents took one of two forms: a Gaussian white-noise process _i, or a short step of current applied at various phases of an ongoing deterministic system. The slowly acting inhibitory conductance g_diD_i was governed by an exponential decay to a cells output firing rate h(V_i) with time constant _i

(2)

A similar schema for self-inhibition was used in Pribe et al. (1997), where the slow self-inhibitory current was construed as an additional neuron providing inhibitory feedback.

The firing frequency of a simulated neuron is a function of the membrane voltage and is determined by a nonlinear sigmoid rate function h(V_i). We chose to use a piecewise polynomial function with the following properties (Jung et al. 1996): outputs rates range from 0 to 1; there exists a true threshold (V_i = 0) below which output rate is equal to 0; the function is sufficiently smooth to facilitate currently unpublished bifurcation analyses

(3)

Analysis of network behavior

Our analysis focused on changes in the durations of parts of the step cycle. All three of the networks analyzed (Fig. 1) included two excitatory interneurons that serve as the model’s output, one to flexor motor neurons (E_f) and one to extensor motor neurons (E_e). The fractionation of the step cycle into parts was determined by time intervals during which these neurons were above a burst threshold. To avoid small noise fluctuations near threshold being registered as separate events, burst threshold was set at 0.1 simulated units of voltage. In addition, periods of threshold crossing <0.1 simulated seconds were not considered to be bursts and were ignored.

The flexor burst (F) and the extensor burst (E) were defined as the period of time in which activity in E_f and E_e, respectively, were above the burst threshold (see Fig. 2). A cycle is defined as the time from the beginning of one flexion phase (F) to the beginning of the next flexion phase (F). The two double bursting models have a second period of flexor activity that is cotemporaneous with extension. This splits the extensor period E into three parts. These are E1, which is the length of time from the beginning of the E burst (extensor on alone) until the beginning of the second flexor burst, E2 in which both E_e and E_f are above threshold, and E3 which is the length of time from the end of E2 to the end of extension. Note there may be short periods during transitions when neither neuron is active, so F+E1+E2+E3 is not always equal to the length of a full cycle. In practice, these small periods of time made little difference to the results.

View larger version (31K): [in this window] [in a new window]   FIG. 2. Half-center output (noise level = 0.03). Voltage traces of the 4 interneurons, 2 excitatory (Ee and Ef) and 2 inhibitory (Ie and If), plus slow self-inhibition in the 2 inhibitory neurons (slow Ie and slow If). The flexion burst (F; solid vertical lines) and extension burst (E; dashed vertical lines) are determined by activity in excitatory interneurons Ef and Ee being above threshold (=0.1). Ie suppresses extension and is active during F; If suppresses flexion and is active during E.

The network architectures studied are shown in Fig. 1. Final parameter values were set with two basic aims: the burst durations should match the burst durations measured during fictive walking in the cat, and overall cycle variability should be similar across models (Boothe and Cohen 2002). The need to control burst durations had to be balanced with the need to achieve stable behavior in the face of added noise.

Overall, the half-center model and clock model (Fig. 1, A and B) are more robust than the feedback model (Fig. 1C), giving functional outputs for many parameters varying over an order of magnitude. A number of parameters affect total cycle duration on both models but keep the basic activity patterns intact. The frequency of oscillation can be increased by raising the tonic drive to the inhibitory interneurons, g_ti, and/or by reducing the magnitude of either g_di or the time constant _i of the slow self-inhibition. Burst durations remain at a fixed relative phase during these changes.

In the clock model, the main difficulty in achieving stable output was a tendency for the model to have two simultaneous winners in the central square of inhibitory interneurons. This was counteracted by increasing tonic excitatory drive to these neurons while increasing the strength of the mutually inhibitory connections (see following text). Because parameters affecting a single interneuron tended to change the durations of multiple phases of the cycle, achieving the appropriate fractionation of the step cycle required the simultaneous adjustment of several parameters.

In the feedback model, a mismatch between parameters affecting the override neuron and those affecting slow inhibition in the central inhibitory neurons I_e and I_f caused nonfunctional alternating outputs where flexion and extension overlap at the end of flexion and the beginning of extension. Due to the hierarchical nature of this model, it was relatively easy to independently adjust the durations of different parts of the step cycle.

Values of the parameters used are reported in the following text. The rest conductance of all neurons in all models was set to g_rI, g_rE = 3.5. Synaptic conductances not represented in Fig. 1 and not described in the following text have been set at , including for all simulated neurons.

The parameters of the half-center model were as follows: the tonic excitation to the excitatory interneurons, g_tE = 1.6; tonic excitation to the inhibitory interneurons, g_tI = 1.0; the conductance of the slow self-inhibition, g_dI = 100; the conductance of the synaptic current from the inhibitory interneurons to the excitatory interneurons, ; and the synaptic current between the mutually inhibitory interneurons, ; and the time constant of self-inhibition on the inhibitory interneurons, _I = 10.

The parameters of the clock model were as follows: the tonic excitation to the excitatory interneurons, g_tE = 2.0; the tonic excitation to the inhibitory interneurons, g_tI = 20; the conductance of the slow self-inhibition, g_dI = 75; the conductance of the synaptic current from the inhibitory interneurons to the excitatory interneurons ; and the medium synaptic current between the mutually inhibitory interneurons, g_syn^I1I4, g_syn^I2I1, g_syn^I3I2, ; the strong synaptic current between the mutually inhibitory interneurons, g_I1I3^syn, g_syn^I2I4, g_syn^I3I1, ; and the weak synaptic current between the mutually inhibitory interneurons, g_syn^I1I2, g_syn^I2I3, g_syn^I3I4, ; and the time constants of self-inhibition, _I1 = 2, _I2 = 3, _I3 = 3, and _I4 = 4.

The parameters of the feedback model were as follows: the tonic excitation to the excitatory interneurons, g_tE = 1.6; the tonic excitation to the inhibitory interneurons, g_tI = 3.5; the tonic excitation to the override neuron, g_tO = 0.7; the conductance of the slow self-inhibition on the mutually inhibitory interneurons I_e and I_f, g_dI = 100; the slow self-inhibition on the override neuron, g_dO = 150; the conductance of the synaptic current from the inhibitory interneurons to the excitatory interneurons, g_syn^IE = 100; the synaptic current between the mutually inhibitory interneurons, ; the synaptic current between the flexor excitatory interneuron and the override neuron, ; and the synaptic current from the override neuron (I_o) to the mutually inhibitory interneurons, , the of self-inhibition on the mutually inhibitory interneurons, _I = 10, _O = 9.

Simulations

Stochastic simulations were run on a Dec-alpha using xpp.aut (Bard Ermentrout, www.math.pitt.edu/bard/xpp/xpp.html). Differential equations used in the model were solved numerically using the Euler-Maruyama method for solving stochastic differential equations (Oksendal 2000). Simulations were originally run using a variety of time steps and methods with little change in model outputs; reported simulations used a time step dt = 0.004. Noise takes the simplest possible form, an additive current into the membrane voltage. At each time step of width dt, noise values in the term _i were drawn from a Gaussian distribution with SD = 1 and then scaled by the quantity /(dt). The final value for was set so that the variability in cycle length was approximately the same for all models ( = 0.03 for the half-center, = 0.04 for the clock model, and = 0.01 for the feedback model).

The deterministic microstimulation experiments were performed using the standard ODE solver available in matlab [ode23; Mathworks, Natick, MA]. Parameters of this model were identical to those used in the stochastic case. A small current was injected into each model for 10 time steps or 0.04 simulated seconds. Due to their greater stability, the magnitude of injected current was larger for the half-center and clock model (0.3) than for the feedback model (0.1). For each of the central inhibitory interneurons in each model, we performed 60 separate simulations in which microstimulation was initiated at regular 0.1-s intervals spanning the 6-s-long prototypical cycle (where s is dimensionless). Burst durations were then compared with perturbation-free simulations. Because perturbations of burst durations were small during the second cycle after the perturbation was given, we analyzed changes in burst duration for phases of the step cycle containing the perturbation as well as the phases in the next cycle.

Measuring and predicting correlations

Correlations between the lengths of different phases of the step cycle were quantified using Pearson’s correlation coefficient, r. For any two variables X and Y, the correlation coefficient ranges from –1 to +1 and is defined by

where Cov(X,Y) is the covariance of X and Y and Var(X) = Cov(X,X) is the variance of X.

Correlation coefficients measured during the simulations with noise could be predicted from the results of the microstimulation simulations described in the preceding text. We write the magnitude of burst duration perturbations as _B (t, I). The subscript B refers to the burst phase for which perturbations are measured (F and E for the half-center model and F, E1, E2, and E3 for the double bursting models), the argument t refers to the time of the perturbation relative to the beginning of the cycle containing the burst, and the argument I is the label of the interneuron to which the stimulation was delivered (I_e and I_f for the half-center; I₁, I₂, I₃, and I₄ for the clock model; I_e, I_f, and I_o for the feedback model). Because all perturbations decay to very small values by the second cycle after stimulation, we only computed values for times in the same and previous cycle as the burst in question, for all other times t we assumed _B (t, I) = 0.

We view noise as a continuous bombardment of perturbations having amplitude a(t,I) delivered at time t to interneuron I. To the degree that perturbations are linear, this perturbation should change the lengths of bursts X and Y by amounts proportional to a(t,I) t_X (t,I) and a(t,I) t_Y (t,I), respectively. Implicit in this linearity assumption is the fact that negative perturbation have an equal and opposite effect as positive perturbations of the same magnitude. This assumption was examined for a sample of perturbations and found to hold to a reasonable degree.

We can perform a prediction pCov(X,Y) for the noise-driven covariance of burst phases X and Y by noting that a given perturbation will contribute to the covariance of X and Y in proportion to the product of the effects of that perturbation on X and Y. Suppose X and Y lie in the same cycle and consider a perturbation delivered to interneuron I at time t relative to the start of that cycle. Let a(t,I) be the amplitude of the perturbation. This perturbation will contribute to the covariance of X and Y in proportion to a²(t,I) t_X (t,I) t_Y (t,I). A prediction pCov(X,Y) can be obtained by averaging these products over the distribution of noise amplitudes and summing these products over all perturbations that affect the lengths of both bursts. Because the noise injections are the same for all interneurons and uniform in time, the average of a²(t,I) is a fixed constant A. Then pCov(X,Y) is proportional to _Itt_x(t,I)t_y(t,I). If burst Y is in the cycle following burst X, then the expression becomes

where the time t is expressed relative to the cycle containing X and L is the length of an entire cycle. Finally

Note that in calculating the values for pCov(X,Y), we sum over time at the resolution at which we delivered stimulation, i.e., time t covers the cycle in 0.1-s intervals. To the degree that this resolution is sufficient to capture the behavior of the system, sampling time more finely should change the preceding expression by a constant factor. Because both this constant factor and the noise scaling A enter equally in the numerator and denominator, they do not affect our final predicted value for the correlation coefficient r_pred.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
We have analyzed three models of walking behavior. We begin with a detailed analysis of the simple half-center model to clarify and highlight basic behavioral principles that emerged from our analysis. We then explore two conceptually distinct mechanisms for the creation of double bursting each of which gives different predictions of the correlation structure between phases of the step cycle.

Half-center model

We start by exploring the behavior of a simple half-center oscillator because it forms the basis for understanding fundamental properties of the double-bursting networks presented later. The half-center model has been studied previously in the context of models of lamprey swimming (Bem et al. 2003; Buchanan 1992; Williams 1992). Dynamics of the half-center are characterized by the combined interplay of the fast winner-take-all interaction between the two mutually inhibitory interneurons (denoted I_e and I_f in Fig. 1A), and the slow dynamics of self-inhibition. Build up of slow inhibition in the currently active neuron eventually tips the competitive balance from one neuron to the other, resulting in alternating activity in the two mutually inhibitory interneurons (Fig. 2). The output of this central CPG is read out by corresponding excitatory interneurons (E_e and E_f in Fig. 1A). Note that because the actions of the central interneurons I_e and I_f are inhibitory, I_e is active during F (flexion) and I_f is active during E (extension).

To examine the dynamic structure within the network, we ran simulations in the presence of a white-noise current injected into each of the two main inhibitory neurons I_e and I_f (see METHODS) and measured overall variability and the correlation between the lengths of extension and flexion bursts in thesame and surrounding two cycles (Fig. 3). Correlation was quantified using the Pearson’s correlation coefficient (see METHODS). This correlation analysis reveals that F is positively correlated with adjacent parts of the step cycle (significantly different from 0; error bars show 2 SD) and is uncorrelated with parts of the cycle further away. This means that a lengthening of F, rather than being correlated with a shortening of the adjacent E’s, actually is correlated with their lengthening as well.

View larger version (12K): [in this window] [in a new window]   FIG. 3. Correlation structure for the half-center model (noise level = 0.03). , Pearson correlation coefficients between flexion burst duration and other bursts for 2 cycles before and after. Error bars on correlation coefficients denote 2 SD. Subscript denotes distance in cycles (0 = current cycle, +1 = next, –1 = previous). F0 indicates the correlation of F with itself (coefficient = 1). E–1 > 0 and E0 > 0 indicate that the flexion burst is positively correlated with the length of the previous and subsequent extension. , linear prediction from microstimulation experiments.

To understand the dynamic mechanisms underlying the correlation structure in the noisy case, we performed a series of microstimulation experiments in which single small perturbations were made to an otherwise noise free simulation (0.2 amplitude current for 0.04 simulated seconds). Figure 4A shows the reaction of the model to microstimulation of interneuron I_e delivered starting at time 2.5. The main effect of the microstimulation is to elongate the flexion burst F. This experiment was repeated, with microstimulation delivered at times evenly spaced throughout the step cycle. The changes in the flexion burst length (F) were then plotted as a function of the time at which the stimulation of I_e was delivered (Fig. 4B). Thus we see that perturbation of I_e late during F results in a relatively large lengthening of that phase, whereas stimulation early during F results in a moderate shortening.

View larger version (18K): [in this window] [in a new window]   FIG. 4. Effects of microstimulation of Ie on the duration of F. a: voltage of Ef in the unstimulated condition (solid line) and with microstimulation starting at time 2.5 (dashed line). Vertical lines mark the end of F in the unstimulated (solid) and stimulated (dashed) simulations. The change in burst duration F = 0.3 s. b: F as a function of the start time for microstimulation. Each point (0.1-s intervals) shows F for a single simulation, e.g., stimulation shown in a contributes the point at 2.5 s. Stimulation of Ie early in the phase shortens F, whereas late stimulation lengthens F.

View larger version (18K): [in this window] [in a new window]   FIG. 5. Effects of microstimulation of Ie (A) and If (B) on burst durations for the current and next cycle. a: voltage traces of excitatory interneurons Ef (solid) and Ee (dashed). Vertical lines show limits of F (solid) and E (dashed). b and c: change in F (F) and E (E) for the cycle undergoing microstimulation. A, a and b are the same as Fig. 4 a and b. Stimulation of Ie near the end of F lengthens both F (time = 1.1 A, b) and E (time = 1.1 in A, c) in the current cycle. Stimulation of If near the end of F shortens F (time = 1.1 in B, b) and E (time = 1.1 B, c) in the current cycle. d and e: change in burst durations of F and E in the cycle following microstimulation (labeled next F and next E above). Stimulation of Ie during E shortens E in the current cycle (time = 2.1 in A, c) and F in the next cycle (time = 2.1 in A, d). Stimulation of If during E lenthens E in the current cycle (time = 2.1 in B, c) and F in the next cycle (time = 2.1 in B, d). Influence of microstimulation on next E is negligible (A, e and B, e).

The effects of microstimulation on the half-center oscillator can be understood by considering three basic properties of networks using fast mutual inhibition and slow self-inhibition: the dynamics are most sensitive to perturbation near the fast competitive events that happen at transition points in the step cycle; the role for a given interneuron switches at these transition points; and microstimulation of a given interneuron has a relatively large immediate influence in one direction followed by a rather long-lasting inhibition that acts in the opposite direction.

The switch from excitation to inhibition is most clearly seen when examining the duration of the burst in which the microstimulation is applied. If the active interneuron is stimulated near the end of a burst, the fast resulting excitation extends that burst (stimulation of I_e extends F, Fig. 5A, b at time 1.2; stimulation of I_f extends E, B, c at time 2.2). This effect is relative large. If the stimulation of the active neuron comes near the beginning of the burst, the direct excitatory effects decay away and self-inhibition hastens the transition to the next phase of the cycle, shortening the ongoing burst (Fig. 5A, b at time 0.6; B, c at time 1.6).

Application of the second property suggests that these effects should be reversed for stimulation of the opposite neuron. This is indeed the case for stimulation late in the burst phase: the fast excitation of the inactive interneuron favors this neuron in the upcoming competition, shortening the burst (Fig. 5A, c at time 2.2; B, b at time 1.2). However, for stimulation occurring early in the burst phase, microstimulation of the inactive neuron does not push this neuron above threshold. As a result, no additional self-inhibition is recruited and stimulation has no effect on burst duration (Fig. 5A, c at time 2.1; B, d at time 1.1).

An important property of the slow self-inhibition is that a relatively brief stimulation can have effects that linger for multiple phases of the cycle. For stimulation late in the burst, the increase in activity that helps the neuron win the competition in the current burst leads to self-inhibition that hinders the neuron during the next cycle. For example, stimulation of I_e near the end of F lengthens that phase but also leads to greater inhibition within I_e. This increased inhibition hinders I_e during the subsequent competition and consequently lengthens E (Fig. 5A, b and c, at time 1.1). Note that in this case the reversal from fast excitation to slow inhibition is accompanied by a reversal in role from the active to the inactive neuron during the transition from F to E. Thus stimulation of the active neuron near the end of a burst increases the duration of both the ongoing and subsequent burst. The same combination of a reversal from excitation to inhibition and a role reversal happens for stimulation of the inactive neuron near the end of the burst, so this stimulation acts to shorten both the current and subsequent burst (Fig. 5A, c and d, at time 2.1 and B, b and c, at time 1.1).

Finally, we point out that the time course over which self-inhibition exerts its influence can be influenced by dynamic events that are influenced by the initial stimulation. For example, consider stimulation of the active neuron during the middle of a burst phase. The stimulation recruits additional self-inhibition that ends the ongoing phase early. However, this reduces the period during which self-inhibition accumulates, largely compensating for the initial increase in inhibition. As a result, after the competitive transition that ends the phase the system has returned to near its unperturbed trajectory and the effect of the perturbation on subsequent phases is negligible (Fig. 5A, b–e, at time 0.1 and B, c–e, at time 1.7).

Correlation structure of the half-center model

Taken as a whole, the microstimulation plots shown in Fig. 5 can be used to understand burst length correlations in the fully noise-driven model (Fig. 3). For example, consider the effects of microstimulating I_e near the end of F. As explained in the preceding text, this stimulation will serve to elongate both F and E. That is, this common cause will tend to induce a positive correlation in the lengths of F and E. Viewing continuous noise perturbation as the sum of many small perturbations distributed in time, we can predict the net correlation between any two parts of the step cycle by computing the product of the effects of microstimulation on each part of the cycle, and summing these effects across all microstimulations (see METHODS). Such a prediction assumes that the effects of perturbing different neurons sum linearly, as do the effects of perturbing the same interneuron at different times. It also assumes that effects scale linearly with the size of the perturbation, and that injecting a negative current has the opposite effect of injecting a positive current. Figure 3 demonstrates that departures from these linearity assumptions are relatively minor.

Clock oscillator model

Like the half-center, the clock oscillator is constructed from a group of mutually inhibitory interneurons acting in a winner-take-all fashion. In generalizing the alternating half-center to a four phase output, the number of inhibitory neurons is increased from two to four (labeled I₁ to I₄, Fig. 1B). The weights within the network are biased such that activity patterns proceed in a single direction around the clock, I₁ to I₂ to I₃ to I₄ (Fig. 6). Each part of the cycle corresponds to activity within a single interneuron. In this sense, the central set of interneurons acts like a clock, marking off the four parts of the cycle. As in the half-center model, there are two excitatory interneurons E_f and E_e providing output to the flexor and extensor motor neurons. Activation of neuron I₁ corresponds to part F (flexor on alone), activation of I₂ to part E1 (extensor on alone), activation of I₃ to part E2 (double bursting), and activation of I₄ to E3 (extension alone).

View larger version (28K): [in this window] [in a new window]   FIG. 6. Clock model output (noise level = .04). a and b: vertical lines show limits of F and E2 (solid) and of E1 and E3 (dashed) in the output of the excitatory interneurons (Ef and Ee). c–f: activity in the 4 central inhibitory interneurons I1, I2, I3, and I4 correspond to phases F, E1, E2, and E3 in serial order.

As a first step in analyzing output variability in the clock model, we set the noise level so that the coefficient of variation (CV = SD divided by the mean) in overall cycle length was 1.5% (noise level = 0.04 gave CV of the cycle = 1.45%). We then measured the CV for each part of the step cycle. All parts of the cycle showed similar degrees of variability as measured by the CV (F: 3.48%; E1: 2.50%; E2: 3.65%; E3: 2.78%). If the parameters of all neurons and connections were exactly symmetric, then the CVs would be identical (up to sampling noise).

We also measured the correlation for each part of the step cycle with all parts of the step cycle one cycle before and one after (Fig. 7). Each of the four subplots in Fig. 7 shows the correlation with a given phase of the step cycle. The figure shows (---) the correlations predicted from applying the linearity assumption to the perturbation results presented in the following text. The basic pattern of correlation is similar for all four phases of the step cycle, i.e., the curve in each of the subplots has roughly the same shape when centered over the phase in question. This means that the strength of the correlation between any two phases of the step cycle is determined mostly by the number of intervening phases of the cycle rather than on the specific phases being compared. Again this is to be expected from the approximate rotational symmetry of the underlying clock network. Second, the strongest correlations are between phases of the cycle that are separated by two intervening phases. This effect is moderately strong, yielding correlation coefficients in the range of 0.05–0.2 (Fig. 7).

View larger version (16K): [in this window] [in a new window]   FIG. 7. Correlation structure for the clock model with noise shown as and —; error bars indicate 2 SD. Each plot shows correlation of 1 phase of the cycle [F (a), E1 (b), E2 (c), E3 (d)] with other phases in the current, next, and previous cycles (subscript = 0, +1, –1, respectively). All phases show the same basic pattern of correlation. Parts of the cycle separated by 2 intervening bursts have the strongest correlation, e.g., F0 is most strongly correlated with E30, E10 with the next F+1, etc. and ---, linear predictions from microstimulation.

View larger version (31K): [in this window] [in a new window]   FIG. 8. Effect of microstimulation of I2 (A, a–i) and I3 (B, a–i) in the clock model. Stimulation I2 and I3 lead to similar patterns of effects relative to the phase in which neuron is active [E1 (A, c) and E2 (B, d), respectively]. Effects are similar to the half-center model. Stimulation of the active interneuron late in its active phase causes lengthening of that phase (time = 2.8 in A, c and time = 4 in B, d); early stimulation causes shortening (time = 2.4 in A, c and time = 3.5 in B, d). Late stimulation also delays the onset time of activity in this interneuron in the next cycle which causes a lengthening of the previous phase (time = 2.8 in A, f and time = 4 in B, g).

The main difference between the models is that the clock model contains a greater variety of functional role reversals at competitive transition points between bursts. In particular, a neuron that has just been active does not participate in the competitions that determine the onsets of the next two burst phases. Therefore the measurable output of the oscillator is not affected until the next cycle, where self-inhibition delays the onset of the corresponding burst. The delayed onset results in a lengthening of the burst before the neuron becomes active. The link between the extension of a burst and the delayed onset of the same burst in the next cycle results in a positive correlation between the length of that burst and the length of the preceding burst in the next cycle (Fig. 7). Negative correlations between remaining bursts are a consequence of mutual inhibition between the inhibitory interneurons of the network.

Feedback oscillator

The feedback oscillator is also a generalization of the basic half-center circuit. But rather than incorporate more neurons within a single winner take all dynamics, this model retains a single pair of mutually inhibitory neurons that control flexor and extensor motor neurons. The basic alternation between flexion (F) and the period of extension (E1-E3) results from the mutual inhibitory dynamics between these two neurons. Double bursting is caused by a separate inhibitory neuron called the override neuron, that inhibits both I_f and I_e (I_o; Fig. 1C). High activity in this cell shuts down activity in all central inhibitory neurons. This releases all excitatory interneurons from inhibition, leading to co-activation of flexor and extensor motor neurons. The override neuron receives excitatory input from E_e. When E_e becomes active at the beginning of extension, it excites I_o, which in turns inhibits I_f leading to a second flexion burst (I_e is already inhibited). This marks the end of E1 and the beginning of E2. Self-inhibition then accumulates in I_o, releasing I_f from inhibition, ending the second flexion burst and terminating phase E2. Voltage traces from the feedback model are shown in Fig. 9.

View larger version (30K): [in this window] [in a new window]   FIG. 9. Feedback model output (noise level = 0.01). Vertical lines show limits of F and E2 (solid; a–e) and of E1 and E3 (dashed; a–e). Top 2 traces show excitatory [Ef (a) and Ee (b)] voltages. Bottom 3 traces show activity in inhibitory interneurons If (c), Ie (d), and Io (e). Io becomes active at the beginning of extension, creating the 2nd flexor burst (time = 2.2 in e). Note the small period of Ie activity near the onset of E3 (time = 3 in d). This is a 3rd period of competition between Ie and If, but one where Ie always loses. This adds significant complexity to the dynamics of the feedback model.

Figure 10 shows the pattern of correlation between different phases of the step cycle. As in Fig. 7, each subplot shows the pattern of covariance with a particular phase of the cycle. In contrast to the clock oscillator, each phase of the step cycle shows a distinct pattern of covariance. Of particular note is the correlation between E3 and the next (F₊₁), which reaches almost 0.6 (Fig. 10, top and bottom).

View larger version (15K): [in this window] [in a new window]   FIG. 10. Correlation structure of the feedback model shown in and with — with noise error bars indicating 2 SD. Phases are labeled as in Figs. 3 and 7. Each individual phase of the step cycle has a specific pattern of correlation with the other parts of the step cycle. Note the very strong correlation between F0 and E3–1, and the different patterns of correlations between the graphs from top to bottom. and ---, linear predictions from microstimulation.

We describe the behavior of the feedback model in the context of how the model responds to microstimulation (Fig. 11), with a focus on each phase of the cycle, starting with the flexion phase F. The transition from F to E1 is governed by the competitive dynamics between I_e and I_f and hence is influenced by the same mechanisms important in the half-center model. In particular, stimulation of I_e (the active neuron) late in the phase extends F, whereas stimulation early recruits self-inhibition and shortens F (Fig. 11B, b). Conversely stimulation of I_f (the inactive neuron) near the end of F shortens F (Fig. 11A, b). The only departure from the behavior of the half-center is that stimulation of I_e near the end of F does not result in a lengthening of the subsequent phase of extension, but rather has a slight shortening effect on E3. This results from a complex interaction between slow self-inhibition in I_e and other dynamic events occurring during extension and will be explained below. Note that the override neuron I_o is below threshold during F, so small stimulations of this neuron have no effect.

View larger version (31K): [in this window] [in a new window]   FIG. 11. Effects of microstimulation on the feedback model. Stimulation of the central inhibitory interneurons If (A, a–i) and Ie (B, a–i) has similar effects on extension and flexion as in the half-center model (compare with Fig. 5). Stimulation of I0 (C, a–i) has the largest effect on the duration of the double-burst phase E2 (C, d) but also has an effect on the transition from E3 to F (C, e). Note that stimulation of all 3 neurons leads to a similar pattern of deviations in the lengths of E3 and the next F (A, d and e, B, d and e, and C, d and e), accounting for the strong positive correlation between these phases (Fig. 10, top, column E3–1, and bottom, column F+1).

The double-burst phase, E2, is terminated by the increasing self-inhibition in I_o that reduces inhibition to the central half-oscillator, allowing I_f to escape from inhibition and end E2. Stimulation of I_o shows the familiar pattern. Stimulation near the end of E2 extends activity in I_o and hence extends E2 (Fig. 11C, d at time 4.5), whereas stimulation of I_o early during E2 recruits greater self-inhibition in I_o and shortens E2 (Fig. 11C, d at time 3.7). If we look at I_f, the effects are in the opposite direction and significantly smaller (Fig. 11A, d). A close examination of the voltage traces during the transition from E2 to E3 (Fig. 9, d) reveals that I_e also increases its activity at this time, initiating another bout of competition between I_e and I_f and causing a dip in the output of E_e. I_f always has the upper hand in this competition (Fig. 9, c), and activity in E_e quickly recovers to support the rest of phase E3. However, the time it takes for I_f to escape from suppression by the override neuron is influenced by competition with I_e. As a result, microstimulation results in a typical half-center pattern, with stimulation of I_e late in the burst extending E2 and stimulation early shortening E2 (Fig. 11B, d). Conversely, stimulation of I_f late in the burst shortens E2 while stimulation early has little effect (Fig. 11A, d).

Finally, we turn our attention to the other transition governed by competition between I_e and I_f: the transition between E3 and F. Again we see the familiar half-center pattern near the end of E3 with stimulation of I_f (the active neuron) extending E3 and stimulation of I_e (the inactive neuron) shortening E3 (Fig. 11, A, e, and B, e). We also see that stimulating I_o near the end of E3 extends E3. I_o inhibits both I_f and I_e (Fig. 11C, e), and so one would not necessarily expect I_o to bias the competition between them. However, a small increase in inhibition has a relatively minor impact on the already active neuron (I_f), whereas I_e is just beginning its rise in activity and the same small increase in inhibition has a relatively large impact. So although I_o sends increased inhibition to both I_f and I_e, it has a larger impact on I_e, delaying the onset of F. If we look further back in time, we see that the length of E3 is also highly sensitive to stimulation of I_e and I_f near the end of E2 (Fig. 11, A, e, and B, e). Interactions between all three interneurons contribute to ending both of the phases E2 and E3. The sensitivity of the dynamics at these two points in the cycle explains the relatively large effects of microstimulation at these times, as well as the fact that phase E3 is relatively more variable than the other phases in the cycle.

As pointed out in the preceding text, the length of E3 is also reduced by stimulation of I_e near the end of F (Fig. 11B, e). Note that this shortening runs counter to expectations derived from the half-center model: stimulation of I_e near the end of F should recruit more self-inhibition in I_e and thus delay the onset of the next F, thereby lengthening E3. The reversal in the effect of slow inhibition is due to an interaction between the slow inhibition and the complex competitive event at the end of E2. With greater inhibition, I_e is less active during this event and this has a net effect of reducing the self-inhibition and hence shortens E3. This illustrates how interactions between long-lasting inhibition and subsequent dynamic events can lead to complex dependencies that span multiple phases of the step cycle.

The correlations seen in the structure in the noise-driven condition (Fig. 10) are entirely consistent with the microstimulation experiments. In particular, note that perturbation of all three interneurons has a similar effect on the length of E3 and the next F. Given the sensitivity of E3 to perturbation, this results in a very large correlation between these two phases. The microstimulation results also reveal that dynamic influences can work in opposition. For instance, for each of the mutually inhibitory interneurons, perturbation during E2 affects the length of both E2 and E3 in the same direction, leading to a positive correlation. However, perturbation of the override neuron in latter half of E2 lengthens E2 significantly, but acts to shorten E3 (Fig. 11, A–C, d and e). This induces a negative correlation. The fact that the net correlation between E2 and E3 is negative indicates that the influence of I_o on this correlation is stronger than the combined influence of I_e and I_f.

Microstimulation effects on extension length in double-bursting models

Differences between the output of the feedback and clock models can also be seen in the influence of microstimulation on the burst duration of the whole of extension (E; Fig. 12). These results are interesting in light of recent results published by Saltiel and Rossignol (2004a,b) (see DISCUSSION). When comparing changes in the duration of extension (E) across models under the effects of microstimulation, three differences are apparent.

View larger version (15K): [in this window] [in a new window]   FIG. 12. Change in total length of extension (E) during microstimulation of interneurons in both the clock (A) and feedback (B) models. For the clock model (A), stimulation of the interneurons affects E in the same manner, i.e., stimulation of I3 during E and I4 during E always shortens E first and then lengthens it, whereas stimulation of the interneurons of the feedback model affects E in a unique way, i.e., stimulation of I0 first lengthens E during E2 and later lengthens E during E3, whereas stimulation of If produces the opposite effect.

Second, in the clock model, microstimulation affects the total length of E in a symmetric manner, i.e., stimulation of any of the three interneurons active during E leads to a part of the phase where E is initially shortened followed by a phase where E is lengthened (compare Fig. 12A, a to c). The same type of stimulation of the feedback model leads to multiple different responses (Fig. 12B). For instance stimulation of interneuron I_o shortens E if delivered at the beginning of E2, elongates E when delivered at the end of E2, and lengthens E when delivered during E3 (Fig. 12B, c). Compare this with microstimulation of I_e where stimulation during E2 acts like stimulation of I_o, but stimulation during E3 gives the opposite effect.

Third, within the feedback model there are times of the step cycle when stimulation of any of the interneurons of the model can influence the burst duration of E, (Fig. 12B, a–c, times 3.2–5.1). Within the clock model there are points in the step cycle perturbable by more than one interneuron (Fig. 12A, b and c, at time 2.8) but these make up only a fraction of the entire cycle, and there are no times where all of the three interneurons influence the burst duration of E.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION GRANTS ACKNOWLEDGMENTS REFERENCES

 
The current study explored three models of the spinal CPG for locomotion. Fundamental properties of networks based on fast competitive inhibition and slow self-inhibition were explored in the context of a simple half-center model. For this model, we find that excitatory perturbation of the active interneuron within a burst has differing effects depending on when in the burst it is given: 1) perturbation early in the burst causes a moderate shortening of the ongoing burst and does not affect subsequent burst durations, and 2) perturbation late in the burst causes a relatively large lengthening of the ongoing burst and a moderate delay in the onset of the corresponding burst in the next cycle. We also find that 3) excitatory perturbation of the inactive interneuron is only effective near the end of the burst and serves to shorten both the current burst and subsequent bursts. Effects 2 and 3 cause adjacent burst durations in the half-center model to either lengthen or shorten together, resulting in a positive correlation between these durations in the simulations with noise (Fig. 3). This positive correlation results from an alignment of the reversal between initial excitation and delayed self-inhibition with the competitive reversal experienced by the stimulated interneuron at the boundary between bursts.

We then generalized the half-center, considering two architectures that give double-bursting output in flexor muscles. In the clock model, the two interneurons of the half-center are replaced by a ring of four mutually inhibitory interneurons with activity proceeding around the ring. The main predictions related to this model stem from the symmetry of the underlying architecture. In particular, all phases of the output have similar CVs, and the strength of correlation between any two phases depends chiefly on the number of intervening phases rather than the particular phases being compared. In the feedback model, a delayed feedback loop interacts with the mutual inhibition underlying flexion and extension to produce double-bursting outputs. The main predictions of this model are that the E1 phase should be relatively insensitive to noise (low CV), whereas phase E3 should have relatively high noise sensitivity (high CV) and that correlations between the burst durations will not depend directly on the temporal relationship between the bursts being compared but on which bursts are being examined as well.

Class of models studied

From a strictly theoretical perspective, the current models expand our understanding of how noise affects the output of neuronal networks that act like relaxation oscillators (Bem et al. 2003; Skinner et al. 1994; Traven et al. 1993; Wang and Rinzel 1995). These models are characterized by large state changes at specific times during the cycle and have distinct fast and slow state variables that provide an explicit substrate for dynamic influences spanning multiple burst phases. The other main approach to modeling spinal CPGs has been to look at the class of oscillators with phase-difference coupling (Kiemel and Cohen 1998; Kiemel et al. 2003). Such phase oscillators are characterized by a single variable (the phase) and generally have state changes that are more continuous.

Although the slow self-inhibition in our models is most commonly attributed to neuronal adaptation via voltage- or calcium-dependent potassium channels (Grillner 1999; Rinzel and Ermentrout 1989), slow inhibition could result from synaptic mechanisms such as G-protein-coupled inhibition and/or from synaptic facilitation in mutually inhibitory synapses. Because the exact temporal profile of the self-inhibition is unlikely to qualitatively alter model behavior, we expect the qualitative predictions derived from our models to apply to a broad class of networks where mutual inhibition and slow self-inhibition drive the dynamics (Morris and Lecar 1981).

Phase-dependent and temporally extended responses

We have shown that the correlation between burst lengths in our models can be accurately predicted by integrating results obtained from microstimulating the interneurons that make up the simulated CPG. In the half-center model, the direct effect of stimulating either central inhibitory interneuron matches the most common action of proprioceptive feedback onto the mammalian CPG (Hiebert et al. 1996; Rossignol and Drew 1988; Whelan 1996). For example, stimulating either Ia or Ib receptors of a given muscle during locomotion tends to enhance activity in the stimulated and active muscle and inhibit activity in muscles antagonistic to them. Stimulation of the inactive muscle (e.g., flexor Ia receptors during extension) can end the current phase early and reset the step cycle to begin the next phase (Hiebert et al. 1996; Whelan 1996). In our model, these results follow directly from the half-center architecture in which flexor and extensor inhibitory interneurons compete for activation.

In addition to the immediate effect on the interneurons driving the oscillation, microstimulation of our models has extended effects that are mediated by slow self-inhibition. For example, the build-up of inhibition in our half-center model leads to the curious property that microstimulation of the active interneuron early in a burst results in a shortening of the current burst duration. A similar difference between early and late perturbation has recently been reported for phasic retraction of a cat’s shoulder during fictive locomotion (Saltiel and Rossignol 2004b; see also Duysens 1977; Duysens and Stein 1978). A similar effect that spans adjacent phases of the step cycle can be seen within the stumble corrective reflex in the context of stumble prevention (Forssberg 1979; Quevedo et al. 2005a). Transient stimulation to the superior peritoneal nerve during extension serves to elongate the extension ipsilaterally and to induce an increase in burst duration of the subsequent flexion (Quevedo et al. 2005a). Our half-center model shows a similar increase in extension and the following flexion after stimulation of I_f near the end of extension (Fig. 5A, b and c). The fact that both phases are elongated is due to the temporal alignment between the reversal from excitation to slow inhibition in I_f and the change in I_f ’s state from active to inactive at the onset of extension.

Because long-lasting effects in our half-center model are sustained by slow self-inhibition, short- and long-term effects must be functionally opposed. Thus our model cannot explain long-lasting alterations in the CPG that consistently favor extension or flexion. For example, our models cannot explain the fact that stimulation of the sural nerve during flexion elongates that flexion burst and also shortens the following extension phase (Duysens 1977) or that stimulation of the superior peritoneal nerve during flexion elongates the current flexion but has no effect on the length of the following extension (Quevedo et al. 2005b).

Transition points in the dynamics of the double bursting models

In a pair of recent papers, Saltiel and Rossignol (2004a,b) argue for the existence of discrete "critical points" within the locomotor output that correspond to the dynamics within bursts rather than the edges of bursts alone. These critical points occur at times when an animal must coordinate changes in activity within multiple muscle groups and are associated with changes in the way the CPG responds to phasic perturbations. Although outputs in our model are restricted to a single pair of antagonistic muscles, the underlying dynamics have two properties that are similar to critical point behavior. First, the major transitions in the dynamics result from strong interactions between interneurons controlling different outputs. As a result, it is possible to alter when these transitions occur, but it is impossible to decouple them without changing the basic behavior of the models. Second, like the results of Saltiel and Rossignol (2004a,b), these coupled events often demark changes in how the models respond to microstimulation. Figure 12 shows the effects of perturbations on the length of extension in our double-bursting