## Abstract

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

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

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).

### Neuron model

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

_{di}D

_{i}, where

*D*ranges between 0 and 1. Voltage dynamics are based on the following differential equation (1) The terms on the right hand side of

_{i}*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*

_{di}D

_{i}and reversal potential

*V*= −1; and a perturbation current

*I*. Perturbation currents took one of two forms: a Gaussian white-noise process ξ

_{i}_{i}, or a short step of current applied at various phases of an ongoing deterministic system. The slowly acting inhibitory conductance

*g*

_{di}D

_{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.

### Model parameters

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. 1*C*), 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 d*t* = 0.004. Noise takes the simplest possible form, an additive current into the membrane voltage. At each time step of width d*t*, noise values in the term ξ_{i} were drawn from a Gaussian distribution with SD = 1 and then scaled by the quantity σ/√(d*t*). 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_{1}, I_{2}, I_{3}, and I_{4} 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 *p*Cov(*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*^{2}(*t,I*) Δ*t _{X}* (

*t,I*) Δ

*t*(

_{Y}*t,I*). A prediction

*p*Cov(

*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*

^{2}(

*t,I*) is a fixed constant A. Then

*p*Cov(

*X,Y*) is proportional to ∑

_{I}∑

_{t}Δt

_{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

*p*Cov(

*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

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. 1*A*), 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. 1*A*). 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.

### Microstimulation of the half-center model

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 4*A* 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. 4*B*). 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.

This technique is generalized to examine how microstimulation of a given interneuron during one part of the step cycle affects the lengths of that and all subsequent parts of the cycle. Figure 5*A*, b–e, shows the effects of microstimulation of I_{e} on the length of four different burst phases: F and E during the current cycle and F and E during the next cycle (Fig. 5*A*, b is an extended version of Fig. 4*b*). Figure 5*B*, b–e, shows the analogous results for stimulation of I_{f}. Notice that the plots in Fig. 5, *A*, b, and *B*, b, fall to zero for times after the end of F simply because microstimulation delivered after the end of F cannot affect its length. A similar causal relationship can also be seen at the right-hand side of Fig. 5, *A*, c, and *B*, c: microstimulation delivered after the end of E cannot affect the length of E. Note also that the symmetry of the underlying model is reflected in the perturbation results. For example, perturbations of I_{f} have the same effect on the length of E as do perturbations of I_{e} on the length of F.

The curves shown in Fig. 5 are related to the well-known phase response curve (PRC) commonly used to examine the behavior of coupled oscillators (Rinzel and Ermentrout 1989; Winfree 2001). However, the PRC measures the effect of a perturbation on the timing of the entire oscillation, whereas Fig. 5 separates out the effect of a given perturbation on the separate subparts of the oscillation over multiple cycles.

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. 5*A*, 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. 5*A*, 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. 5*A*, 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. 5*A*, 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. 5*A*, 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. 5*A*, 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. 5*A*, 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_{1} to I_{4}, Fig. 1*B*). The weights within the network are biased such that activity patterns proceed in a single direction around the clock, I_{1} to I_{2} to I_{3} to I_{4} (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_{1} corresponds to part F (flexor on alone), activation of I_{2} to part E1 (extensor on alone), activation of I_{3} to part E2 (double bursting), and activation of I_{4} to E3 (extension alone).

Dynamics within the clock model have the property of rotational symmetry. For example, interneuron I_{1} has a similar relationship with its neighbors I_{4} and I_{2} as interneuron I_{2} has with its neighbors I_{1} and I_{3}. More generally, rotating the labels in Fig. 1*B* would not qualitatively change the pattern of relationships between the four central interneurons of the model. Because activity in each interneuron is associated with one burst phase of the cycle, we expect that the relationship of phase E1 with the surrounding phases F and E2 should resemble the relationship of phase E2 with its surrounding phases E1 and E3, etc. Note that because the time constants for the slow self-inhibition were altered to roughly match the experimentally measured burst lengths, the underlying symmetry of the dynamics is only approximate.

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).

We explore the underlying dynamics of the clock model by examining the effect of microstimulation on the central inhibitory neurons driving the oscillation. Figure 8 shows the effect of microstimulation on interneurons I_{2} and I_{3}. These plots again reveal the underlying rotational symmetry of the clock model, with stimulation of neuron I_{2} having the same qualitative pattern as stimulation of I_{3}, just shifted in time and phase.

All the main features of the microstimulation experiments can be understood from direct application of the three properties outlined for the half-center model in the preceding text. Stimulation of the active neuron results in the same large lengthening when delivered near the end of the burst and a moderate shortening when delivered earlier (Fig. 8, *A*, c at times 2–2.8, and *B*, d at times 3–4). Like the half-center, the shortening of the phase from early stimulation largely compensates for an increase in self-inhibition so that early stimulation has little effect beyond the current phase. Also, stimulation of an inactive neuron near the end of a given phase aids this neuron in the upcoming competition and hence shortens that phase (Fig. 8, *A*, b at time 1.9, and *B*, c at time 2.9).

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. 1*C*). 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.

Unlike the clock model, the feedback model is hierarchical, and different mechanisms contribute to the transition between different phases of the step cycle. In particular, the transitions from F to E1 and from E3 to F are determined by mutual inhibition between I_{e} and I_{f}, the transition from E1 to E2 is determined by the buildup of activity in I_{o}, and the transition from E2 to E3 is determined by the decay of inhibition from I_{o} as well as the ability of I_{f} to escape from this inhibition. As a result of this heterogeneity, the different parts of the step cycle show very different degrees of variability. Again we adjusted noise level so that the CV of cycle length approximated 1.5% (noise level σ = 0.01 gave CV of the cycle = 1.71%). Unlike the clock model, the different parts of the step cycle in the feedback model have widely different CVs (F: 2.97%; E1: 1.09%; E2: 3.19%; E3: 4.86%).

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_{+1}), which reaches almost 0.6 (Fig. 10, *top* and *bottom*).

### Microstimulation in the feedback oscillator

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. 11*B*, b). Conversely stimulation of I_{f} (the inactive neuron) near the end of F shortens F (Fig. 11*A*, 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.

Next we examine the initial phase of extension, E1. Looking at the voltages in Fig. 9, we see that the onset of extension is marked by a rapid rise of activity in the excitatory extensor interneuron E_{e} (Fig. 9, b, time 2), followed with a short delay the rapid rise in the override neuron I_{o} (Fig. 9, *e*). This causes a rapid decline in I_{f} activity (Fig. 9, *c*), leading to double bursting. Note that all transitions are strong and rapid. That means that even though perturbations may have a significant impact on voltage, these changes will be converted into rather minor differences in the time at which threshold is crossed. Consistent with this, stimulation of all interneurons, including the override neuron I_{o}, has a very small effect on the length of E1, explaining why E1 is the least variable phase in this model.

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. 11*C*, d at time 4.5), whereas stimulation of I_{o} early during E2 recruits greater self-inhibition in I_{o} and shortens E2 (Fig. 11*C*, d at time 3.7). If we look at I_{f}, the effects are in the opposite direction and significantly smaller (Fig. 11*A*, 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. 11*B*, d). Conversely, stimulation of I_{f} late in the burst shortens E2 while stimulation early has little effect (Fig. 11*A*, 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. 11*C*, 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. 11*B*, 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.

First, microstimulation of some interneurons in the feedback model affect the length of E across a large portion of E (microstimulation of I_{o} in Fig. 12*B*, c), whereas the same stimulation in the clock model only alters E if delivered in a small portion of the burst (microstimulation of I_{3} in Fig. 12*A*, b).

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. 12*A*, a to c). The same type of stimulation of the feedback model leads to multiple different responses (Fig. 12*B*). 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. 12*B*, 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. 12*B*, 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. 12*A*, 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

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. 5*A*, 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 models. The length of extension is most sensitive to perturbation near the dynamic transitions associated with the double-bursting events.

### Coordination of complex movements

The complexity and time course of the effects within our models suggest how transient inputs could aid in the coordination of complex movements. For example, transient sensory input to the spinal cord might trigger a series of dynamic changes that alter muscle activity over time periods extending into the next step cycle. Alternatively, temporally localized inputs from higher centers could trigger chains of functionally related effects to coordinate voluntary movements such as obstacle avoidance (Drew et al. 2004). Such long-lasting effects from transient stimulation have been observed under conditions of very strong stimuli from motor cortex (see Fig. 10 in Bretzner and Drew 2005) and may explain why a complicated behavioral repertoire remains in animals without extraspinal inputs (e.g., Forsberg 1979). Of course not all coordination will be accomplished within the spinal CPG, as it is well known that both motor cortex and sensory information modulate the CPG on a cycle-by-cycle basis (reviewed in Cohen and Boothe 2002; Drew et al. 2004).

Such extended changes in burst lengths due to perturbation are often attributed to the triggering of complex activity patterns within excitatory and inhibitory interneurons (Bretzner and Drew 2005; Quevedo 2005b). Our simplified models raise the possibility that some of these functional reversals may reflect the shift from transient excitation to the dominance of slow self-inhibition. Moreover, our double-bursting models demonstrate that simple switches from immediate excitation to long-lasting inhibition can have complex and surprising effects that depend on the functional architecture of the model. In the clock model, for example, each interneuron only participates in the competitions at the boundaries of a single phase of the cycle. Therefore perturbation of a single interneuron changes the length of the current burst and has very little effect on the two intervening bursts and then another substantial effect on burst length. Somewhat surprisingly, this temporally separated interaction is not seen between the same phase across cycles, but rather a given phase in one cycle and the preceding phase in the next cycle. This is because burst durations are not coded per se but rather are determined by competition-driven state transitions at the beginning and end of each phase. A different form of complexity is demonstrated by the feedback model, in which a short period of activity in one of the central interneurons (I_{e}) at the end of phase E2, has a very minor effect on the output of the CPG but serves to ‘reverse the sign’ of the effect of slow inhibition on later parts of the step cycle (see results).

### Clock-like versus feedback model of double bursting

Much of the research on spinal CPGs addresses two basic questions. What controls the timing relations between various motor outputs? Is the same circuit or different circuits used to accomplish qualitatively different behaviors? The phenomenon of double bursting allows one to study a blend of these questions within the same basic behavior. In particular, it raises the question of whether the same or different parts of the spinal CPG drive the two bursts of flexion. A more specific but related question is whether the CPG is dominated by a clock-like mechanism. If the dynamical mechanism is simply marking out time, one expects an underlying symmetry in the dynamics in which distance in time rather than similarity in peripheral effect, is the dominant variable determining the dynamic structure of the network (cf. Golubitsky et al. 1999).

Although quantitative predictions cannot be drawn from the abstract models considered here, the models do suggest a number of qualitative predictions. Most generally, a clock hypothesis predicts that correlation patterns for all phases of the cycle should be qualitatively similar (e.g., Fig. 7). More specifically, in a double-bursting system the phases of flexor activation (F and E2) lie on “opposite sides” of the step cycle. Under a clock hypothesis, these phases are expected to be “maximally different” and have little correlation. One also expects that perturbing the CPG during these two phases should lead to quite different effects. In contrast, in a model in which the two phases of flexor activation are driven by the same part of the CPG, one would expect commonalities between these two phases, such as perturbations having similar effects in phase F and E2 (e.g., compare the effect of stimulating I_{e} and I_{f} during F on ΔF Fig. 11 and stimulation during E2 on ΔE in Fig. 12).

### Measuring correlations as well as perturbations

A classic method for studying a system is to perturb it and observe the pattern of effects. We have shown that noise-driven variability can also reveal useful information about the structure of the underlying system. In our models, stimulating distinct populations of interneurons within the CPG cause perturbations that are generally consistent with these noise-driven correlations. In fact, the noise-driven correlations can be predicted by combining the perturbations from all components that drive the oscillation. In the actual system, however, little is known about the correlation structure of burst lengths during fictive walking and whether particular physiological perturbations act in concert with or in opposition to the natural variability within the system. It also is an open question of whether perturbations, in addition to their effect on mean burst length, act to change patterns of subsequent variability in the system. Our approach suggests that analyzing the trial-to-trial variability under a number of experimental conditions may shed light on the structure of the underlying CPG.

## GRANTS

This work was partially supported by National Institute of Neurological Disorders and Stroke Grant NS-39909 to A. H. Cohen and a Sloan Research Fellowship to T. W. Troyer.

## Acknowledgments

The authors thank Dr. Tim Kiemel for helpful comments.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 by the American Physiological Society