## Abstract

In earlier work, we have developed an integrated model for insect locomotion that includes a central pattern generator (CPG), nonlinear muscles, hexapedal geometry and a representative proprioceptive sensory pathway. Here, we employ phase reduction and averaging theory to replace 264 ordinary differential equations (ODEs), describing bursting neurons in the CPG, their synaptic connections to motoneurons, muscle activation dynamics and sensory neurons, with 24 one-dimensional phase oscillators that describe motoneuronal activation of agonist–antagonist muscle pairs driving the jointed legs. Reflexive feedback is represented by stereotypical spike trains with rates proportional to joint torques, which change phase relationships among the motoneuronal oscillators. Restriction to the horizontal plane, neglect of leg mass and use of Hill-type muscle models yield a biomechanical body–limb system with only three degrees of freedom, and the resulting hybrid dynamical system involves 30 ODEs: reduction by an order of magnitude. We show that this reduced model captures the dynamics of unperturbed gaits and the effects of an impulsive perturbation as accurately as the original one. Moreover, the phase response and coupling functions provide an improved understanding of reflexive feedback mechanisms.

## 1. Introduction

Legged animals display great agility in traversing noisy, uncertain environments, but their neuro-muscular dynamics and control strategies remain poorly understood. They also provide interesting examples of hybrid dynamical systems: piecewise-holonomic constraints (Ruina 1998), due to intermittent leg contacts with the ground, imply that their equations of motion switch as stance patterns change or flight phases occur (Holmes *et al.* 2006). Insects are of particular interest since they are small and fast; their feet contact the substrate so briefly that reflexive neural feedback typically cannot act within a single stride, suggesting that their stability might primarily derive from preflexive mechanisms due to mechanical reaction forces (Brown *et al.* 1995; Holmes *et al.* 2006).

Here, we develop a simplified hybrid dynamical model of insect locomotion in the horizontal plane, using observations on the death’s head cockroach *Blaberus discoidalis* (Full & Tu 1990, 1991; Full *et al.* 1991; Jindrich & Full 2002), and drawing on an earlier lateral-leg-spring model (Schmitt & Holmes 2000*a*,*b*; Schmitt *et al.* 2002) that represented each of an insect’s stance leg tripods with a single passive spring, analogous to the spring-loaded inverted-pendulum description of sagittal plane dynamics (Blickhan 1989; Ghigliazza *et al.* 2003). Subsequently, Seipel *et al.* (2004) introduced six actuated legs, which were later improved to include joints actuated by torsional springs (Kukillaya & Holmes 2007) and then by muscles (Kukillaya & Holmes 2009). Ghigliazza & Holmes (2004*b*) modelled the central pattern generator (CPG) and motoneurons, and Kukillaya *et al.* (2009) assembled an integrated neuro-mechanical model, driven by a CPG and equipped with rudimentary reflexive feedback.

Such models have shown that feed-forward actuation, coupled with piecewise-holonomic constraints, can provide partial asymptotic stability (Ruina 1998; Holmes *et al.* 2006), consistent with the preflexive stabilization hypothesis, prompting experimental studies that confirm this (Jindrich & Full 2002). Animals are nonetheless endowed with many proprioceptive and exteroceptive sensors, and we wish to use our models to reveal how these reflexive control circuits interact with the feed-forward dynamics.

Kukillaya *et al.* (2009) considered feedback of leg forces to motoneurons, and Proctor & Holmes (2008) studied simple steering strategies, showing that these could enhance resistance to, and recovery from, perturbations. However, the model of Kukillaya *et al.* (2009), containing some 270 ordinary differential equations (ODEs) describing bursting neurons and synapses, the calcium release cascade of muscle activation and sensory neurons, is too complex to permit direct analyses. Here, we use phase response curves (PRCs) and averaging theory to reduce the key motoneuronal dynamics to a small set of ODEs that capture phase relationships and describe their modulation by tonic feedback of joint torques.

The paper is organized as follows. In §2, we review the integrated neuro-mechanical model of Kukillaya *et al.* (2009), further details of which can be found in Kukillaya & Holmes (2009), Kukillaya (2010) and the electronic supplementary material. We outline the reduction method and compute PRCs and averaged coupling functions in §3; these describe how the CPG and motoneurons interact and react to force feedback, and hence how muscle activation timing is set and adjusted during running. Section 4 presents the results of simulations and analysis, and a discussion follows in §5. Background on legged locomotion models appears in Holmes *et al.* (2006) and in numerous other works cited below.

## 2. The neuro-mechanical model

Here, we review the muscle-actuated hexapedal model developed in Kukillaya & Holmes (2009) and the CPG–motoneuron model of Ghigliazza & Holmes (2004*b*), and describe their union in the integrated model of Kukillaya *et al.* (2009). Further details can be found in the electronic supplementary material and other references given in this section.

### (a) Biomechanical components: body, limbs and muscles

Many insects, including cockroaches, employ a double-tripod gait over a wide speed range (Full & Tu 1990, 1991). The left front and rear and the right middle legs support the body in the left (L) stance phase; the right front and rear and left middle legs then take their place in right (R) stance, while the former tripod swings forward. Building on the hexapedal models of Seipel *et al.* (2004) and Kukillaya & Holmes (2007), Kukillaya & Holmes (2009) introduced muscle models with nonlinear length and velocity dependence (Zajac 1989), excited by pulses representing motoneuronal spikes. In these and the present model, agonist–antagonist muscle pairs actuate ‘hips’ and ‘knees’ modelling the coxa–trochanter–femur and femur–tibia joints active in horizontal plane motions (body–coxa joints that raise and lower the legs, and minimally actuated tibia–tarsus joints, are excluded). As in earlier work, we assume a 50 per cent duty cycle: each tripod’s touchdown (TD) coincides with liftoff (LO) of the contralateral tripod.

Figure 1 illustrates the model’s geometry in relation to the cockroach. As in the earlier work, a single rigid body models the head, thorax and abdomen, and leg mass (approx. 5–10% of total mass) is neglected, yielding three degrees of freedom: two translational and one rotational. Since the dynamics is invariant with respect to translation in the plane and LO and TD events coincide, the system state at TD of each tripod is specified by four variables: centre-of-mass (CoM) velocity magnitude *V* , velocity direction *δ* relative to body axis, body axis angle *θ* and angular velocity , and the dynamics is determined by an L-TD to L-TD Poincaré map (Seipel *et al.* 2004; Kukillaya & Holmes 2007):
2.1

In Kukillaya & Holmes (2009), appropriately phased pulse inputs were computed to produce body dynamics characteristic of cockroaches by solving an inverse problem to determine the required foot forces and joint torques, and then using an optimization routine to select appropriate pulse timings. The ‘target’ foot forces, detailed in electronic supplementary material, S1, are only used to derive pulse timing.

Driven by these pulses, the body obeys Newtonian equations of motion during each stance phase, written here in inertial (ı,j) ground-plane coordinates (figure 1*b*):
2.2
where denotes the reaction force at foot *i*, the vectors **r**_{fi} identify TD foot positions (assumed fixed during each stance phase) and sums are taken over *i*=1,2,3 for the L tripod and *i*=4,5,6 for the R tripod (see numbering convention of figure 1*b*). The foot forces are obtained from joint torques via nonlinear kinematic relationships inherent in the geometry of figure 1*b*: see electronic supplementary material, S2.

Joint torques include passive (linear) stiffness and damping and moments due to muscle forces, described by relations of the form
2.3
(Hill 1938; Zajac 1989), where *F*_{0} is a scale factor and *A*(*t*) an activation function that depends on the input pulses and models the release of calcium and its binding to contractile fibres (Hatze 1977, 1978). (The *F*_{0} are selected in the inverse fitting procedure.) The functions *F*_{L} and *F*_{V} describe force dependence on muscle length *l* and velocity *v* and, hence, on joint angles and angular rates, so mechanical feedback enters here; see electronic supplementary material, S3, and Kukillaya & Holmes (2009) for further details. Equipped with the joint torques, equation (2.2) forms a closed, nonlinear, hybrid dynamical system with prescribed time-dependent forcing and TD/LO events, which determines foot forces and body motions.

Appropriately parametrized, this model possesses a partially asymptotically stable branch of gaits with a single unit eigenvalue corresponding to body orientation *θ* (the heading direction). The other three eigenvalues lie within the unit circle throughout the relevant speed range, and two of them remain small in magnitude, implying that, following a brief transient of approximately one-step duration, the dynamics of recovery from perturbations is essentially two-dimensional, evolving near a subspace spanned by the two larger eigenvalues (Kukillaya & Holmes 2009).

### (b) Neural components: the central pattern generator and motoneurons

Little is known regarding locomotor centres in cockroaches, but the three thoracic ganglia contain bursting interneurons, and it is believed that neighbouring ganglia have ipsilateral inhibitory coupling, and contralateral inhibitory coupling may also exist (Pearson & Iles 1970, 1971, 1973; Pearson 1972). Based on these and more recent studies of stick insects (Büschges 1998; Büschges *et al.* 2004, 2008), a CPG model with six bursters, each driving slow and fast bursting extensor motoneurons in a feed-forward manner, was constructed by Ghigliazza & Holmes (2004*a*,*b*) using single-compartment, ion-channel ODE models (figure 2).

Each of the six CPG units represents a hemi-segmental ganglion involving multiple interneurons. These establish the basic double-tripod pattern and stepping frequency, units 1, 2 and 3 bursting in phase, and 4, 5 and 6 likewise, but 180^{°} out of phase with 1, 2 and 3. Following Pearson & Iles (1973), extensors (active primarily in stance) are periodically inhibited by the CPG units, and so burst approximately out of phase with them, while flexors (active in swing) are excited and burst in phase with the CPG. (The tonically active extensors support the animal when standing still with CPG inactive.) The CPG burster and motoneuron models are briefly described in electronic supplementary material, S4.

### (c) The integrated feed-forward model

In the integrated model of Kukillaya *et al.* (2009), the pulse inputs were replaced by motoneuronal spike trains as triggers for muscle activation. The animal must coordinate many more than the two pairs of muscles per leg of the model, and multiple CPG interneurons may produce detailed timing information during the stepping cycle to coordinate distinct slow and fast muscles actuating each joint. The pulse sequences derived to fit foot forces in Kukillaya & Holmes (2009) are quite varied, but in spite of the fact that the CPG produces only two stereotypical bursts per full L–R stride (cf. equation (3.8)), it was nonetheless possible to tune individual motoneuron bias currents and conductances to approximate those sequences and produce acceptable foot forces. Below, we use the averaged coupling functions to parametrize the reduced system in a more transparent fashion (cf. §3*b* and the figures therein).

### (d) Proprioceptive feedback of joint torques

The exoskeletons of insect legs contain sense organs called campaniform sensilla (CS), whose neurons fire at rates approximately proportional to strain. The locations and varied orientations of CS render them sensitive to force direction (Zill & Moran 1981*a*,*b*; Zill *et al.* 1981), so that magnitudes and signs of joint torques *τ*_{±} can be monitored.

Rather than modelling feedback pathways associated with each joint by sensory neurons with spike rates proportional to positive and negative torques in that joint as in Kukillaya *et al.* (2009) and Proctor & Holmes (2010), here we assume that tonic spike rates are determined directly by torque:
2.4
with baseline rate *b* (Hz) and constant of proportionality *n* (Hz N^{−1} m^{−1}). Feedback is assumed present only during stance while the leg transmits force, and is relayed, via excitatory and inhibitory synapses, to extensor and flexor motoneurons activating the muscles of the joint in question, so that increased resistance to leg extension is compensated by increased torque generation and vice versa (figure 3). Positive torque means that extensor force exceeds flexor force, tending to extend the leg, and vice versa, the directions being opposite in knees and hips (figure 1*b*). (This will be generalized in §4 to allow feedback relative to an average torque level in unperturbed running.)

Feedback from each leg to motoneurons actuating other legs is probably present (Borgmann *et al.* 2009), and sensory cells may also synapse on CPG neurons, thereby influencing the underlying stepping rhythm, but we exclude such complications in the present model. Further details on sensor models and other feedback pathways appear in Proctor & Holmes (2010).

## 3. Phase response curves, averaging and a reduced model

We now outline phase reduction theory (Malkin 1956) and show how it and the averaging theorem (Guckenheimer & Holmes 2002, §§4.1 and 4.2) are used to simplify the CPG and the motoneuron model of figure 2. The theory strictly applies only to weakly perturbed oscillators, but in practice it often gives good indications of behaviour for strongly coupled systems.

### (a) Phase response curves

We assume that the bursting model of §2*b*, detailed in electronic supplementary material, S4, equations (S21)–(S27), possesses an attracting hyperbolic limit cycle *Γ*_{0} (Guckenheimer & Holmes 2002). A local coordinate system can be defined in a neighbourhood of *Γ*_{0}, with a phase variable *ϕ*∈*S*^{1} that follows *Γ*_{0} at a constant frequency *ω*_{0}, absent external perturbations and synaptic inputs. The equations (S21)–(S22) governing each neuron take the form
3.1
where **x**=(*v*,**w**) contains the membrane voltage and gating variables, **f** denotes the intrinsic vector field and *ϵ***g** contains all input currents (except the fixed bias that establishes bursting). The structure of the flow along *Γ*_{0} and exponential attraction towards it implies that the state space is locally foliated by smooth transverse *isochrons* of dimension *n*, such that any two solutions starting in one leaf, , at time *t*_{0} will lie in another, , at time *t*_{1} (Guckenheimer 1975). Solutions starting on any therefore approach *Γ*_{0} with the same asymptotic phase, implying that we may write *ϕ*=*ϕ*(**x**) and **x**=**x**(*ϕ*) and, using the chain rule, deduce that
3.2

The phase variable *ϕ* is normally defined so that *ϕ*=0 marks the start of the burst (e.g. maximum voltage of the first spike). Since input currents enter only via the first (voltage) component *g*_{v} of **g**, equation (3.2) becomes
3.3
where the PRC characterizes the ‘receptivity’ of the system to inputs as a function of phase. Figure 4 shows an example of a PRC for a bursting neuron (cf. Ghigliazza & Holmes 2004*b*, fig. 7).

### (b) Averaged evolution equations for weakly coupled systems

From equation (3.3), each reduced phase equation may be written as
3.4
where *g*_{ji} and *h*_{si} denote coupling from other bursters and inputs from sensors to neuron *i* and the *ϕ*_{i} dependence in these functions is due to post-synaptic voltage dependence (electronic supplementary material, S4). Rapid phase evolution is removed by defining slow phases *ψ*_{i}=*ϕ*_{i}−*ω*_{0}*t*, and each component of the vector field is then averaged over the unperturbed period *T*_{0}=2*π*/*ω*_{0} and the terms dropped:
3.5
The averaging theorem implies that solutions of equation (3.5) remain within 𝒪(*ϵ*) of those of equation (3.4) on a time scale of , and that hyperbolic equilibria of equation (3.5) correspond to hyperbolic periodic orbits of equation (3.4) (Guckenheimer & Holmes 2002, theorem 4.1.1). Here, we assume that the CPG sets the stepping frequency, so that motoneurons are entrained and the time (*t*) dependence in *h*_{si}(*ϕ*_{i},*t*) is *T*_{0}-periodic. Letting *τ*=*ψ*_{i}+*ω*_{0}*t*, the integrals of equation (3.5) become
3.6
revealing that the coupling function *G*_{ji} depends only upon phase differences (*ψ*_{i}−*ψ*_{j}), but that feedback terms *H*_{si} can introduce the absolute phase *ψ*_{i}.

For two identical units symmetrically coupled with no external inputs (*g*_{12}=*g*_{21}=*g* and *h*_{si}≡0), the averaged system takes the form
3.7
and the phase difference *θ*=*ψ*_{1}−*ψ*_{2}=*ϕ*_{1}−*ϕ*_{2} satisfies
3.8
Since *G*(0)−*G*(−0)=0 and *G*(*π*)−*G*(−*π*)=*G*(*π*)−*G*(*π*)=0 (because *G* is 2*π*-periodic), in-phase and anti-phase solutions *always* exist for this symmetric system. Their stability is determined by the slopes and 𝒢′(*π*).

As explained in Ghigliazza & Holmes (2004*b*), the bilateral symmetry of the CPG circuit (figure 2) allows it to be reduced to a pair of phase-locked units and thus to the simple form of equation (3.8), where *θ*=*ϕ*_{L}−*ϕ*_{R} denotes the phase difference between the left and right tripods. The coupling function is found to qualitatively resemble sin*θ*, so that *θ*=*π* is the unique stable fixed point (Ghigliazza & Holmes 2004*b*, fig. 8). Thus, the CPG forms a clock that sets the stepping frequency by emitting two equi-spaced bursts of spikes per cycle to entrain motoneurons that coordinate the left and right leg tripods.

The motoneurons of figure 2 are excited or inhibited by the CPG, so we compute averaged coupling functions *G*_{ji}(*ψ*_{i}−*ψ*_{j}), where *ψ*_{j} denotes the CPG phase and *j*=L or R for the left and right tripods, respectively, with *ψ*_{R}=*ψ*_{L}+*π*. Examples of coupling functions for two of the 12 motoneurons actuating the six joints of one tripod appear in figure 5. The positions of stable fixed points (*G*_{ji}(*θ*)=0,*G*_{ji}(*θ*)′<0) are primarily determined by the frequency difference between the two units (Proctor & Holmes 2010), as the dashed and solid curves illustrate. Specifically, if the isolated CPG (unit 1) has frequency *ω*_{0} and the motoneuron (unit 2) frequency *ω*_{0}+*ϵδ*_{2}, then the phase equation for the latter takes the form
3.9
which yields the averaged system
3.10
where *ψ*_{i}=*ϕ*_{i}−*ω*_{0}*t* and *τ*=*ψ*_{2}+*ω*_{0}*t* (electronic supplementary material, S4). Thus, changes in detuning *δ*_{2} shift *G*_{12} vertically while coupling strength scales its magnitude, so that entrainment ranges increase approximately linearly with coupling strength (Proctor & Holmes 2010, figs 4 and 5). Since inhibitory (respectively excitatory) inputs produce coupling functions similar to negative (respectively positive) sine waves (figure 5), extensor (respectively flexor) motoneurons should have higher (respectively lower) natural frequencies than the CPG to ensure entrainment. Furthermore, the slope *G*_{ji}(*θ*)′ at fixed points determines not only eigenvalues, but also the sensitivity of fixed-point locations to changes in frequency: there is a stability/flexibility trade-off, with strong stability implying that larger parameter changes, or stronger feedback, are required to modify the motoneuron phase. Here and throughout, we set *ϵ*=0.1.

These facts guide parameter selection. In particular, the averaged coupling function in figure 5*b* shows why it was difficult to obtain phase-locked solutions with appropriate spike timing for some flexor motoneurons in Kukillaya *et al.* (2009). Since flexor connections are excitatory, the motoneuron frequency must typically be lower than that of the CPG to obtain entrainment, but to provide the desired phase shift, it cannot be much lower, and the resulting stable fixed point is therefore close to a saddle-node bifurcation that occurs at a slightly lower motoneuron frequency. Hence, the fixed point is quite sensitive to excitatory inputs, because they add positive terms to the right-hand side of equation (3.10), bringing the fixed points closer to collision. The extensor motoneuron can accept larger excitatory inputs without losing entrainment, and both coupling functions can clearly survive large inhibitory inputs, which contribute negative terms. Phase differences *ψ*_{2}−*ψ*_{1} specify motoneuron timing with respect to L and R TDs, which are assumed set by the CPG.

Proprioceptive feedback introduces terms of the form
3.11
into the system (electronic supplementary material, S4), in which the synaptic input *s*(*t*) depends upon the torque in the joint in question, as further described in §4. In the presence of unexpected perturbations, this cannot be predicted throughout the cycle, and so *h*_{s2} cannot be averaged, but must be added to the vector field of equation (3.10) and integrated directly using the substitution *ϕ*_{2}=*ψ*_{2}+*ω*_{0}*t*. In general, the *j*th motoneuron’s phase is determined by
3.12
where *ψ*_{1}=0 for the L tripod, *ψ*_{1}=*π* for the R tripod and , respectively, denote excitatory and inhibitory feedback terms, which are distinguished by their reversal potentials *E*^{post}_{2,syn}.

Following the analysis in Proctor & Holmes (2010, §3.4), *in unperturbed periodic running*, the feedback terms *h*_{s2} can be further approximated by assuming that tonic spike rates of sensory neurons are sufficiently high that each synaptic input *s*(*t*) behaves as if it were almost constant throughout its firing interval *t*∈[*t*^{0},*t*^{1}]. Letting *τ*=*ψ*_{2}+*ω*_{0}*t*, this yields
3.13
where *Δ* is the average value of *s*(*t*) over [*t*^{0},*t*^{1}], which in turn scales linearly with torque as specified in equation (2.4). Examples of this function appear in figure 6.

## 4. Simulations and analysis

The model developed in §3 contains a CPG clock and 24 scalar ODEs of the form (3.12) for the motoneuron phases in place of the 264 ODEs describing the CPG, motoneurons and sensory neurons in Kukillaya *et al.* (2009). Tonic feedback enters each *h*_{sj} via the time-varying function *s*(*t*) in equation (3.11) that advances or retards the motoneuron phase. Since the full dynamics of the bursting neurons are not represented, spikes cannot be added or deleted, nor interspike intervals modified, as in Kukillaya *et al.* (2009). The simplicity of this model evidently comes at a price, which we assess in the following numerical experiments.

In all simulations, the CPG bursting frequency was set at 9.98 Hz. Each muscle is activated according to equation (2.3), using the analytical form of Kukillaya & Holmes (2009), when the phase variable of its motoneuron passes values corresponding to spikes in the integrated model of Kukillaya *et al.* (2009): see electronic supplementary material, S3, for details. Integration of the three-degree-of-freedom system (2.2) with these inputs produces an average speed of 0.275 m s^{−1} in unperturbed running, within the preferred speed range of *B. discoidalis*, and having variations of roughly ±0.025 m s^{−1}, similar to those observed (figure 7).

### (a) Running at preferred speed

In this section, we compare gaits with and without joint torque feedback. Full model details, including parameter values used for the various simulations, are given in the electronic supplementary material.

#### (i) Feed-forward running: gait and eigenvalues

With sensory pathways disabled, the feed-forward model runs stably, producing the gait with net foot forces, CoM moments and periodic dynamics illustrated in figure 7, in comparison with data taken from a freely running animal. Solving for the fixed point (*V* =0.30, *δ*=0.07, *ω*=4.4 and *θ* arbitrary, angles given in radians), and estimating the linearized Poincaré map as described in Kukillaya & Holmes (2007), we confirm that the eigenvalue corresponding to body orientation is unity, as expected due to rotational invariance, and that the remaining three eigenvalues are less than one in magnitude. As in Kukillaya & Holmes (2009), two eigenvalues are very small (𝒪(10^{−3})), and the third is ≈0.54, so that the strongly stable modes decay within a single L–R stride, and there is a two-dimensional slow subspace tangent to the span of the slow and the neutral eigenvectors. The slowly decaying dynamics in this subspace dominates the recovery from perturbations.

Using the muscle model parameters and spike timings of Kukillaya & Holmes (2009), the model runs somewhat faster than the two samples of data shown in figure 7 (top left), but speeds with the range 0.2–0.3 m s^{−1} are commonly observed. The major difference between model results and the data is in the impulsive forces that occur around TD and LO (right column), but as the figure indicates, integration of the ODEs substantially smoothes these and the resulting velocity and yaw variations are very similar to those observed in experiments. The impulses are presumably due to TD and LO occurring when foot forces are relatively large, and alternative protocols for TD and LO that take forces and joint angles into account may reduce their magnitudes. We are currently studying this.

In spite of the significant yawing oscillations and lateral motions, we refer to gaits like that of figure 7 as nominally straight running, because the asymptotic CoM path is bounded by two parallel straight lines (figure 8).

#### (ii) Running with force feedback: gait and eigenvalues

Feedback enters the motoneuron phase equations via synaptic currents, as described in equations (3.11)–(3.12). The spike rate *f*_{s±} of the pre-synaptic sensory neuron depends on the magnitude of joint torque via equation (2.4), in which *τ*_{±} are computed from the variation around the mean torque over an unperturbed stance. Specifically, assuming that , we write *s*(*t*)=*b*_{1}+*n*_{1} |*τ*_{±}|, where *b*_{1} and *n*_{1} are scaled versions of the parameters *b* and *n* in equation (2.4) that are chosen to keep the (dimensionless) synaptic variable *s*(*t*) in an appropriate range (electronic supplementary material, S4). This is done by first examining torque profiles during unperturbed feed-forward running to establish maxima, minima and the average level 〈*τ*〉, and then applying the impulsive perturbation described below, again noting maxima and minima. This determines effective ranges for *τ*_{±}, bounded below by 〈*τ*〉, and above by the magnitudes of the maximum and minimum perturbed torques, respectively. The parameters *b*_{1} and *n*_{1} are then chosen so that *s*(*t*) occupies a range consistent with that observed in simulations of the full spiking neuron model (Proctor & Holmes 2010), consistent with spike rates observed in unperturbed and heavily loaded animals (Noah *et al.* 2004); see electronic supplementary material, S4, for details. Note that only one of *τ*^{+} or *τ*^{−} is active at a given time during stance.

Summarizing, the ODEs describing each motoneuron contain four terms: the original (constant) frequency *ω*_{0}+*ϵδ*_{j}, the CPG coupling function *G*_{1j} and excitatory and inhibitory coupling from the sensors, as in equation (3.12). Addition of joint-torque feedback changes the coordinates of the fixed point to *V* =0.32, *δ*=0.14, *ω*=3.6, and modifies the gait accordingly. The strongly stable eigenvalues remain of , and the weaker at ≈0.5.

### (b) Responses to an impulsive perturbation

We now repeat the perturbation experiment of Jindrich & Full (2002), confirming that the phase-reduced model correctly describes the insect’s response to a lateral force impulse while running in the preferred speed range. We employ the same triangular impulse of duration 4 ms as that used in Kukillaya & Holmes (2007) and Kukillaya *et al.* (2009), which is an idealized version of the rapid impulse generated by the cannon carried by the insect in Jindrich & Full (2002), delivered approximately 5 ms following left TD.

Figure 8*a*,*b* shows the impulses and the resulting lateral CoM velocities recorded in Jindrich & Full (2002), and produced by the model during a 230 ms period including the perturbation; typical unperturbed velocities in nominally straight running are also shown for comparison. CoM paths in the inertial (ground) plane appear in figure 8*c*. Model results are shown for the feed-forward system and for three cases with torque feedback from all joints, one with only excitatory synapses active and two with both excitatory and inhibitory synapses, as in the full circuit of figure 3. The latter two cases have different inhibitory synaptic conductances, weak (dark grey) and stronger (light grey), but both considerably weaker than the average excitatory conductance (see electronic supplementary material, S4, for values).

In all cases the lateral dynamics of the response are similar to those of a cockroach running at a preferred speed, as shown in figure 8*a*. Lateral velocity magnitudes are comparable, although the rise time is faster for the model and its decay rate is slower. The resulting CoM path and body orientations are also similar to those of the mechanical model, showing rapid preflexive recovery to straight running, but with substantial changes in heading. (Magnitudes of these depend upon the timing of the impulse with respect to TD: cf. Kukillaya & Holmes (2007, fig. 14); data not shown here.) Notably, final changes in heading are significantly smaller for the cases with feedback, and smallest for the case with full feedback and weaker inhibitory conductances.

These results are consistent with the feed-forward simulations of Kukillaya & Holmes (2007, figs 14 and 15) and of Kukillaya *et al.* (2009, figs 8–10), in which only feedback from knee torques was included. They amplify the latter in showing that modulation of muscle forces and consequent changes in active joint torques can partially compensate for the effects of perturbations, and that the degree of compensation depends upon a balance of excitatory and inhibitory effects. In particular, it is interesting that weak inhibition provides greater compensation than strong inhibition, suggesting that there might be an optimal inhibition/excitation ratio.

## 5. Conclusion and discussion

We have developed a simplified neuro-mechanical model of a running insect with tonic proprioceptive feedback of joint torques in legs. The detailed dynamics of the CPG and bursting motoneurons are replaced by scalar equations that preserve key dynamical phenomena while reducing complexity by an order of magnitude. The reduced model yields explicit approximations of CPG–motoneuron phase relationships, thereby permitting analyses of the relative effects of feed-forward and feedback pathways. We show that the model captures key features of the overall system dynamics, including gait characteristics (figure 7), and responses to impulsive perturbations (figure 8).

The analysis of §3 shows that fixed points of phase-locked motoneurons in the feed-forward system depend on two factors: CPG–motoneuron frequency difference and the strength of CPG–motoneuron synapses. As the coupling function’s slope grows with synaptic strength, eigenvalues increase and a greater frequency difference is required to produce a given phase difference. However, while stability is clearly desirable, feedback should be capable of modifying phase relationships in response to sensory inputs, and there is clearly a trade-off between these effects, as described previously (Proctor & Holmes 2010). Weaker stability and a flexible neural system may provide greater resilience to perturbations than strongly stable feed-forward architectures that leave little room for reflexive modification. The fact that the CPG and sensory input terms are additive in the reduced motoneuron model (equation (3.12)) makes predictions of phase changes in steady running due to different levels of feedback rather simple (e.g. consider linear superpositions of figures 5 and 6 for different torque levels and durations). It also allows interpretation of the effects of time-varying feedback in situations such as the impulse experiment of §4 (e.g. consider linear superposition of phase-shifted curves from figure 6 due to positive and negative torques).

Using this model, we have shown that tonic feedback can modify CPG–motoneuron phase relationships and thereby partially counter impulsive perturbations, leading to smaller heading changes following recovery to nominally straight running. We lack the space here to describe systematic studies of the effects of feedback strength, or the balance of inhibition and excitation noted in §4*b*, but this relatively simple and analytically tractable model will enable future studies of this type. We are also interested in relaxing coincident TD and LO of left and right tripods (50% duty cycle), and considering LO rules for individual legs based on joint angles or foot forces as well as the CPG phase (as at present). We note that gaits of the hexapedal models of Kukillaya & Holmes (2007, 2009) have already proved robust to random variations in foot TD positions, TD timing and motoneuron spike timing, so we expect the present model to accommodate the variability that such rules would introduce.

## Acknowledgements

This work was partially supported by NSF EF-0425878 (Frontiers in Biological Research) and the J. Insley Blair Pyne Fund of Princeton University.

## Footnotes

One contribution of 12 to a Theme Issue ‘Theory of hybrid dynamical systems and its applications to biological and medical systems’.

- © 2010 The Royal Society