Insects running on elastic surfaces

Running animals diverse in leg number and posture often negotiate complex, heterogeneous environments (Dickinson et al., 2000). These environments may have complex spatial structure, be composed of materials with non-linear mechanical properties, be dynamic, or more likely show some combination of these attributes. To sustain rapid locomotion, legs must cyclically contact the environment, modulating forces experienced by the body such that task level goals, such as continuous forward progress relative to an external feature, are attained (Koditschek et al., 2004). An animal may, or may not, be able to maintain stability as it attempts to move through a complex, unpredictable environment at a given speed. Cockroaches, for example, maintain forward velocity when confronted with rough terrain three times their hip height (Sponberg and Full, 2008), whereas lizards have been found to slow down, change kinematics and pause when confronted with obstacles (Kohlsdorf and Biewener, 2006). Failure to maintain stability increases the probability of extreme yawing, pitching or rolling, and may result in a reduction in forward speed, all of which are likely to be highly detrimental to the survival of an animal that depends upon rapid running.

General models of running are now being used to explain stability. One class of proposed models consists of those referred to as ‘spring–mass systems’. Driven by the finding that across a broad range of size and leg number, the centre of mass motion of runners qualitatively resembles that of simple, bouncing spring–mass systems (Blickhan and Full, 1993; Farley et al., 1993), these models collapse the anatomical details of individual legs and the action of multiple legs during simultaneous contact to a single, virtual, elastic spring leg. The spring-loaded inverted pendulum or SLIP model describes sagittal plane motion with a point mass atop a linear spring (Blickhan, 1989), and an analogous model of horizontal plane motions of sprawled posture runners is the lateral leg spring or LLS model (Schmitt and Holmes, 2000a; Schmitt and Holmes, 2000b). These models, and their underlying parameters, such as leg stiffness, landing angle and average forward velocity, have been variously used to describe and predict aspects of legged running and hopping. Human runners, for example, have been found to exhibit constant SLIP model virtual leg stiffness against speed (McMahon and Cheng, 1990) and gravity (He et al., 1991), which they achieve mainly through increased landing angle in the former case and reduced vertical landing velocity in the latter. Extending to non-humans, Farley and colleagues examined seven running and hopping animals of various sizes, and found that leg stiffness increased with size, but because landing angle remained constant, larger animals had longer overall contact duration, and thus lower stride frequency (Farley et al., 1993). A functional explanation for the generality of spring–mass-like behaviour is debated; it may supplant enhanced locomotor economy through storage and return of elastic strain energy in structures such as tendons (Alexander, 1988), but the fact that spring–mass mechanics are maintained even when a lossy surface makes it more costly to do so (Lejeune et al., 1998; Moritz and Farley, 2003) suggests that another role, such as stability, can also be important (Geyer et al., 2005).

Studies of running on compliant surfaces, largely conducted on humans, have shown that humans can adapt their leg function when confronted with these surfaces. In some of these studies the aforementioned spring–mass models were used to interpret these changes. McMahon and Greene observed changes in step length and contact time on surfaces of varying stiffness that agreed with a model of the runner as a spring-damper (dashpot) bouncing into a single linear elastic spring model of the surface (McMahon and Greene, 1979). They found that a specific surface stiffness could tune the runner–surface system for faster running speed, through longer stride lengths that outweigh a slightly lower step frequency (McMahon and Greene, 1979) (but see Stafilidis and Arampatzis, 2007). This was followed by the finding that humans can maintain their centre of mass (COM; for a system of particles or a continuous mass distribution the COM is the average of the positions of each piece of mass weighted by the mass at each position; the system's mass can be considered concentrated at this point for many calculations) (see Thornton and Marion, 2004) motion when running on elastic surfaces, which has been interpreted as an adjustment of their SLIP leg stiffness to compensate for the surface (Ferris et al., 1999; Ferris et al., 1998). Kerdok and colleagues explored the mechanism by which this was accomplished, finding it to be a more extended knee posture that causes the limb to drive further into a softer surface (Kerdok et al., 2002). Concomitantly they found that metabolic energy consumption was reduced on the elastic surface, demonstrating that human runners were able to shift some portion of the burden of bouncing the body forward to the surface. When challenged with very soft surfaces, humans have been found to reverse the phasing of their leg movements, extending their leg and driving into the surface upon contact, and flexing the leg as it unloads after mid-stance (Moritz and Farley, 2005). Finally, on lossy and damped surfaces, human runners and hoppers also maintain their COM motion, but the leg transitions to acting as a power-producing actuator (Lejeune et al., 1998; Moritz and Farley, 2003). These results show that human runners and hoppers can maintain COM motion across an array of compliant surface conditions.

A running animal may adjust the neural commands that it sends to its musculoskeletal system to maintain stability when it is confronted with the aforementioned perturbations, or it may continue to drive the system with unaltered commands. The former strategy is referred to as a feedback control strategy, whereas the latter is a feedforward control strategy. An extensive literature exists on feedback control strategies used by stick insects and other arthropods performing slow, quasi-static walking behaviour (Büschges et al., 2008; Cruse et al., 2007; Ritzmann and Büschges, 2007), including the finding that feedback-mediated control strategies change with substrate compliance (Cruse et al., 2004). Feedback is used to control individual legs (Burrows, 1992), limb loads (Noah et al., 2004), and inter-leg coordination of the timing of leg movements (Cruse, 1990). Rapid running behaviour may present a serious challenge to these feedback control strategies, however, because not enough time is available to process and react to perturbations (Jindrich and Full, 2002; Kubow and Full, 1999). A neural feedback control strategy is limited by the time required for sensory transduction, afferent transmission, computation, efferent transmission and muscular force development. As a result it has been hypothesized that context-dependent control is critical to robust locomotion, and that slower behaviours will take advantage of the available time and utilize feedback, whereas faster behaviours will resort to feedforward control (Koditschek et al., 2004).

Successful feedforward control relies on a system that can continue to move stably in the face of perturbations even when driven with control signals that do not respond to these perturbations (Sponberg and Full, 2008). One way to achieve this is for the system to either be or emergently act as a mechanical system that can self-stabilize in the face of perturbations (Kubow and Full, 1999). A self-stabilizing system either absorbs de-stabilizing mechanical energy changes caused by the perturbation or transduces them away from de-stabilizing motions into other modes (Brown and Loeb, 2000; Schmitt et al., 2002). The aforementioned SLIP and LLS models have been shown to be self-stabilizing (Geyer et al., 2005; Ghigliazza et al., 2005; Schmitt and Holmes, 2000a). Self-stabilization means that within certain ranges of their parameter space and for certain types and magnitudes of perturbation, these models can continue to move successfully after being confronted with the perturbation. In addition to relying on the properties of the model system to recover from perturbations, control strategies may be layered on top of the model system, by adjusting parameters of the model during locomotion. Running humans confronted with a variable height, rigid step upwards exhibit reduced virtual leg stiffness on the perturbed step, which can be interpreted as moving within the parameter space of the SLIP model, potentially to a region which is more self-stabilizing (Grimmer et al., 2008).

Self-stabilization often refers to the system as a whole; for constituents of the system the term mechanical feedback has been used to refer to musculoskeletal structures that stabilize locomotion in the face of perturbations without the use of neural feedback. The response of these structures may depend on their own current state or that of the animal, and hence is a form of feedback, but they are passive in the sense that their stabilizing behaviour does not require active, neural control. For example, spines on the leg of the spider Hololena adnexa latch onto debris during leg extension, but collapse passively when the leg is pulled free (Spagna et al., 2007). Thus the purely passive, mechanical performance of the structure is dependent on its state and it provides a form of state-dependent mechanical feedback.

The fact that stabilization is an integrated phenomenon that spans from individual appendages to COM motion raises the question of whether control itself is centralized or decentralized (Koditschek et al., 2004). Intuitively, a centralized control scheme uses a single, central clock signal to drive many appendages, whereas a decentralized control scheme might have separate independently running clock signals for each appendage, that influence each other to a greater or lesser degree. Moving between these schemes is performed by changing the degree to which decentralized control signals are tied together. The importance of the centrally generated pattern has recently been demonstrated by the research of Li and colleagues, who found that small changes in the phase and frequency of the central clock signal have a dramatic effect on the ability of a six-legged robot to run successfully over sand (Li et al., 2009). In bipedal runners confronted with an unexpected, rigid drop-step perturbation, for instance, a proximo-distal gradient in the degree of control exerted has been found (Daley and Biewener, 2006), suggestive of at least some degree of decentralized control. Which organization is most effective will depend on the requirements of the locomotor task, available sensors, their bandwidth and their noise levels, the nature of external perturbations, and the ability of distributed mechanical feedback to handle local perturbations. Research has found that for environments with rigid, but spatially discontinuous supports, many-legged, sprawled posture runners use such distributed mechanical feedback and the bridging effect of kinetic energy to run stably (Spagna et al., 2007; Sponberg and Full, 2008).

The aim of this study was to determine whether multi-legged runners can adjust the stiffness of their virtual leg. If they can do so, we would have important evidence to support the notion that a centralized controller can tune the emergent behaviour of the system through the action of multiple legs. If they cannot do so, then we would seek to understand whether the observed mechanics are the result of an unaltered control strategy. To test this hypothesis, we required a multi-legged runner that exhibits SLIP-like running mechanics, that has documented virtual leg properties, and for which the neural control architecture has been studied. We therefore selected cockroaches, Blaberus discoidalis, because they exhibit SLIP-like running mechanics (Full and Tu, 1990), the stiffness of their virtual leg spring on rigid surfaces has been measured (15 N m⁻¹) (Full and Tu, 1990), and their neural control has been studied during slow (Mu and Ritzmann, 2005; Ridgel and Ritzmann, 2005; Ridgel et al., 2007; Ritzmann and Buschges, 2007; Tryba and Ritzmann, 2000a; Tryba and Ritzmann, 2000b; Watson and Ritzmann, 1998a) and fast (Sponberg and Full, 2008; Watson and Ritzmann, 1998b) locomotion.

We chose an elastic surface perturbation because it presents the simplest possible mechanical element with which to perturb the virtual elastic spring leg, for which an adjustment of the virtual leg can be readily interpreted, and because the response of single stance leg runners to this perturbation is known. We chose a surface stiffness of 10 N m⁻¹, lower than the animal's virtual leg stiffness because the predicted displacement of surfaces stiffer than the leg spring (~2 mm) represents a small perturbation when considering the recently discovered ability of the animal to handle rough, albeit rigid, terrain varying in height by 30 mm (Sponberg and Full, 2008). We define the ability to maintain forward speed on an elastic surface as the existence of running trials in which the mean forward speed on the elastic surface is not significantly slower than that on the rigid surfaces before and after it.

We consider this study as demonstrative of a new horizon in animal locomotor studies, in which microinstrumentation, judiciously and appropriately applied, can avail the investigator of a new class of data (Byrnes et al., 2008). To measure the COM mechanics of the cockroach as it moves over the elastic surface we fabricated a novel, miniature sensor that gives six degree of freedom information about the animal's COM dynamics. The sensor combines microelectromechanical systems (MEMs) inertial sensing of linear acceleration with sophisticated Kalman filter-based automatic video tracking and state estimation. The result is an ability to measure the COM dynamics of intact, freely behaving animals as they move in complex environments, removing the constraint of small area, rigid, planar force plates. The device has been successfully integrated with standard electromyography (EMG) techniques (Sponberg et al., 2007), and future engineering research will combine it with flexible multielectrodes capable of resolving several muscles simultaneously (Spence et al., 2007). The reduced size of this instrumentation is critical to future comparative biomechanical studies of locomotion, as it enables the study of smaller animals with drastically different form that, because of their size, move in very different physical regimes (Spence, 2009). Understanding how they do so will accelerate our discovery of the general principles that govern the neuromechanical control of locomotion.

We used adult Death's head cockroaches (Blaberus discoidalis L.) of both sexes (Carolina Biological Supply, Gladstone, OR, USA). Their average mass was 2.33±0.48 g (mean ± s.d.; N=6 animals). Animals were housed in large plastic containers and fed dried dog food and water ad libitum. Experiments were performed at room temperature (24°C).

Cockroaches with accelerometer backpacks ran across the long axis of a 48 cm×28 cm rectangular Plexiglas arena where the ‘floor’ of the arena was Plexiglas apart from a central rectangular section of dimensions 21 cm×11 cm that was an elastic latex membrane (Fig. 1A). A high-speed video camera (Kodak Ektapro HG2000, Rochester, NY, USA) directly above the arena recorded each trial. On selected trials, a second high-speed video camera recorded a simultaneous side view.

Accelerometer backpacks were designed such that the top surface of the accelerometer package could be attached directly to the cuticle, ensuring that the sensor was placed as close as possible to the COM (Fig. 1B). The electrical signal connector and video tracking cross were mounted on the opposite face of a printed circuit board (PCB), and thus projected dorsally from the animal. Animals were anaesthetized in CO₂ for 2 min, after which the forewings (tegmina) and hindwings were trimmed to expose a square 7 mm×7 mm area of the dorsal cuticle just caudal to the thoracic–first abdominal segment joint. Care was taken not to cut larger circulatory pathways of the wings. The exposed cuticle and accelerometer surfaces were very gently abraded, and the backpack was attached with cyanoacrylate adhesive (LOC30379 Super Glue Gel, Loctite, Rocky Hill, CT, USA), and held in place for 30 s while the adhesive cured. Cockroaches were then allowed at least 1 h to recover at room temperature prior to running trials.

After recovery, the cockroach was placed at the centre of one end of the long axis of the arena, and the electrical tether connected to the backpack. Rapid running behaviour was elicited by gently probing the posterior abdominal segments and cerci with a small rod. Cockroaches ran quickly into the arena, over compliant surface, and encountered the far wall of the arena. If the cockroach stopped, either in the middle of the arena or upon making contact with a side wall of the arena, we repositioned it manually at the centre of the start end of the arena. We defined a successful trial as one in which the cockroach made a complete traversal of the arena, running continuously from the starting point to the opposite wall. After a successful trial, the cockroach had at least 5 min to recover as the video was downloaded from the camera buffer. We continued recording until at least 10 trials meeting this operational definition were obtained.

Each trial was divided into constituent steps, using the phase estimation technique described below. As we wished to test whether cockroaches could maintain forward speed on the compliant surface, some trials, and some steps, were not included in the analysis. We rejected steps under two conditions: (1) when the cockroach turned more than 15 deg., and (2) when it changed speed by more than 0.15 m s⁻¹. We then identified acceptable trials for analysis. To control for history effects using our A–B–A experimental design, we rejected trials in which the step variable of interest was significantly different on the rigid surfaces before and after the compliant surface using a Wilcoxon rank sum test (P<0.05 resulted in rejection). Within the trials meeting this criterion, we considered only those in which at least four acceptable steps (meeting the turn and speed change criteria above) occurred, in sequence, on each of the three surface conditions (rigid surface before compliant surface, compliant surface, and rigid surface after compliant surface). To carry out a balanced statistical analysis with equal numbers of steps per surface type, we compared the sequence of four steps on the rigid surface before the compliant surface with the sequence of steps on the compliant surface. This resulted in a total of 288 steps being analysed: four steps per surface type, two surface types, six trials per animal and six animals.

Forward speed was computed as the mean forward speed across all time points in the step. Step frequency was computed as the reciprocal of the duration, in seconds, of the step. The peak-to-peak dorsoventral acceleration was computed as the difference between the maximum positive and minimum negative peaks of the dorsoventral acceleration during the step. To avoid pseudo-replication in the peak-to-peak dorsoventral acceleration data, in the form of neighbouring steps using the same negative peak as their minimum, the negative peak value was taken from the data occurring within the step, but before the maximum peak only.

Normality of the data was tested with the Shapiro–Wilk test, and accepted data were analysed using a three-factor repeated measures ANOVA. The three within-subjects factors were: step number (1–4), surface (rigid before compliant substrate, compliant substrate), and trial number (1–6). Statistical analyses were carried out in Matlab (The Mathworks, Inc., Natick, MA, USA), using custom-written scripts, the Statistics Toolbox, the RMAOV33 routine (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=9638) and in SPSS (SPSS Inc., Chicago, IL, USA).

We fabricated a compliant substrate by clamping a stretched latex membrane (Sheer Glyde Dental Dams Model No: GLD-200, Glyde USA, Seattle, WA, USA) across the open rectangle in the Plexiglas arena floor. A CNC milled stainless steel press ring was slotted into a matching groove in the Plexiglas floor, uniformly stretching the membrane. A force lever (model 300B; Cambridge Technology, Inc., Cambridge, MA, USA) in force control mode was used to measure the force–length relationship of the membrane by applying linear ramps in force over time (Fig. 1C). Utilized membranes had a stiffness between 8.5 and 13 N m⁻¹ (11.2±2.0 N m⁻¹, mean ± s.d., N=4), and displayed linear force–length relationships (Fig. 1C). Repeated measures of the membrane stiffness over periods of up to 5 days demonstrated that membrane properties were stable for the duration of the experiments. To ensure that the membrane presents a linear elastic substrate to the animal at its step frequency of 25 Hz, and that the acoustic effect of added mass and damping due to deflection of air is negligible, we conducted frequency sweep experiments on membranes made of this material. These sweeps showed a resonance peak at approximately 100 Hz, with a width of 25 Hz, significantly higher than the step frequency of the animal, demonstrating that the membrane will appear linearly elastic to the animal.

We designed a backpack so that it provided a 3-axis linear COM acceleration measurement with a dynamic range of ±2 g, allowed an unencumbered, freely behaving animal to traverse an arbitrary terrain, and had suitable means of validation against ground truth measurements. We met these aims by fabricating miniature backpacks around a 3-axis MEMs inertial sensor (MMA7260, Freescale Semiconductor, Austin, TX, USA), incorporating a five-point retro-reflective marker cross for estimating rigid body dynamics from videography, and interfacing the backpack with a detachable, lightweight, flexible wire tether (Fig. 1B).

The surface mount accelerometer chip was soldered to a small printed circuit board (ExpressPCB, Santa Barbara, CA, USA; www.expresspcb.com), whose other components consisted of a power supply stabilization capacitor (0.1 μF, SMT0603 package, DigiKey Inc., Thief River Falls, MN, USA) and a five-pin male header (Mill-Max P/N 850-10-050-10-001000, DigiKey Inc.). Ground, +3.3 V power, and analogue X, Y and Z acceleration voltages were traced to the header pins. Five conductor micro-tethers were constructed to interface the backpack on the animal with data acquisition electronics. Tethers were fabricated by soldering five 2 ft (0.61 m) long strands of 0.002 in (0.05 mm) diameter epoxy insulated silver EMG wire (California Fine Wire Co., Grover Beach, CA, USA) between five-pin female headers (Mill-Max 851-93-050-10-001000, DigiKey Inc.), and twisting the resultant cable to form a tightly braided tether. The computer side of the tether was connected to an interface PCB that provided the power rail and passive RC filtering (cut-off frequency=1600 Hz) of the accelerometer voltages. Filtered signals were acquired directly into Matlab at a 1 kHz sampling frequency (PCI-MIO-64E-1 board, National Instruments, Austin, TX, USA).

The video tracking cross was constructed of lightweight balsa wood, and consisted of four markers in the plane of the accelerometer and PCB, placed a distance R=1 cm from the centre of the PCB, projecting outward in each of the rostral, caudal, left and right lateral body axis directions. A vertical beam established a fifth, top marker, 3.6 cm above the centre point of the lower cross. This configuration was designed such that a single 2D image from above yields six degree of freedom information about the position and orientation of the cross, and hence the backpack.

Individual accelerometers and the separate axes of each accelerometer had varied voltage offsets and sensitivities. As such, each accelerometer was calibrated by placing the backpack at the centre of a rotating head and spinning the device through 360 deg. while recording the output voltages. Rotating the device in the Earth's gravitational field in this manner presented the device with a known 2 g signal swing, centred about the zero-g voltage offset. First the accelerometer X and Y axes were calibrated by spinning the device about its Z axis, followed by calibration of the Z axis by spinning the device about its X axis. Devices calibrated in this manner exhibited 2% error in their estimate of the magnitude of g when being rotated in place, equivalent to 0.5 mN force plate resolution for a 2.5 gram cockroach. This is well within the limits expected from the manufacturer.

High-speed video data (512×384 pixel frames) were collected at 500 frames s⁻¹. Simultaneous hardware triggering of video and accelerometer data acquisition produced synchronized data. We collected video data to measure body orientation (pitch, roll and yaw), and to measure the position and velocity of the animal. The world-to-camera coordinate transformation was calibrated using a direct linear transformation (DLT) matrix, which was computed from static images of a vertically staggered lattice of bricks (Lego Systems Inc., Enfield, CT, USA) with retro-reflective markers at precise and known x, y and z offsets (Abdel-Aziz and Karara, 1971; Hatze, 1988; Heikkila and Silven, 1997; Hinrichs and McLean, 1995). Image processing, data acquisition, signal processing and statistical analyses were all carried out using custom-written scripts and graphical user interfaces (GUIs) in Matlab.

Video tracking of animal motion using an unscented Kalman filter We implemented an unscented Kalman filter (UKF) for automatic tracking of video sequences and direct estimation of the cross, and hence the animal's, rigid body configuration (position, velocity and orientation within the world frame). The Kalman filter (KF) is a computational tool widely used in tracking, estimation, sensor fusion and control applications, because of its relative simplicity, optimality and robustness (Julier et al., 2000; Julier and Uhlmann, 2004; Wan and van der Merwe, 2000). The UKF is a more recent evolution of the KF that has been shown to be particularly effective at tracking through non-linear functions, while being straightforward to implement.

The KF is a recursive estimation technique whereby a set of observations taken over a discrete set of time points are used to estimate the internal state of a system (Maybeck, 1979). The operator must supply a model of how the system evolves with time, a model of how the system state is mapped to observations, and estimates of the noise present in both the system and observations. In our implementation, the system state is a 10 dimensional vector describing the position and orientation of the cross in a fixed world, or laboratory, coordinate system (the 10 dimensional state vector consists of the orientation angles α, β, γ; the x, y and z position; and the x and y velocity, and x and y acceleration). The z coordinate is assumed to be fixed (the animal is constrained to move on the surface of the arena), and hence we do not include z velocity and acceleration in our system state. We arbitrarily chose the origin of the world x, y and z coordinate system to be the bottom left corner of our calibration object, and angles α, β, γ to specify the orientation of the object in the standard sequence of Euler angles. These angles were converted to pitch, roll and yaw for data analysis and interpretation using trigonometry.

The system evolution function we utilize in the UKF is one of constant acceleration motion; for a time step Δt, x and y position are updated with vΔt + 1/2aΔt², whereas velocities are incremented with aΔt. This does not mean that the filter will not estimate these quantities when passed over observation data; it simply specifies the underlying model the filter uses to predict how the system changes with time. The observation function takes the rotations and translation specified in the state vector, applies them to a 3D model of the cross, and then uses the DLT to compute the corresponding 2D image (x, y pixel) coordinates for each of the five markers of the cross. In essence, the UKF determines the configuration of the cross in the world coordinate state vector that is the most reasonable updated estimate given the previous estimate and the observed 2D marker positions. We utilized the UKF predictions of where the cross markers will be to solve the matching problem in subsequent video frames (i.e. to identify which identified marker in the picture corresponds to which point on the cross). This allowed automatic tracking of the cross through each video, with user input only required in specifying the initial marker locations.

We first validated the accuracy of the DLT camera model by affixing the cross to a brick (Lego) and translating it, with orientation constrained by the brick on a base-plate (Lego), to varied locations spanning the field of view of the camera (Fig. 1D). The UKF tracker was allowed to converge on the static marker positions, and the error in its estimate of the direction of the top cross marker was measured to be 0.57±0.35 deg.

Validation of the accelerometer, video cross tracking and combined signal processing algorithm was carried out using Lego calibration objects, both statically (Fig. 1D) and dynamically, through placement of the calibration object on a rotating turntable. The backpack, affixed to a brick (Lego), was placed on a base-plate (Lego) on the turntable at a precise distance from the spindle and spun at 45 r.p.m. This moved the backpack through a precise, known trajectory in time and space, against which our tracking system was calibrated. The error in measured velocity of the backpack, as compared to the known velocity (computed with ν=ωr; where ν is linear velocity of the backpack, ω is known angular velocity of the turntable and r is radius from the centre at which the backpack was placed, measured using Lego studs) provided by the turntable, was 0.29±1.1 mm s⁻¹ (mean ± s.d.). The error in the system's estimate of gravity at rest was 0.004±0.015 g, and the error in the system's estimate of dynamic acceleration of the backpack (after subtraction of g), against the known centripetal acceleration of a=ν²/r, was −0.020±0.021 g. To estimate how much error arose from sources of DC bias such as thermal drift and power supply irregularity, we fitted the theoretically expected acceleration trajectory to the data obtained on the turntable, and found the error about the theoretical trajectory to be −0.005±0.012 g.

We computed forward and lateral speed by projecting the velocity vector of the COM onto the forward and left lateral cross axes, respectively, and low-pass filtering the result (zero-phase, third-order Butterworth, 100 Hz, used for all quantities to follow except where noted). We calculated pitch as the angle between the z=0 plane (the arena floor) and the forward cross axis (positive upward), and roll as the angle between the floor plane and the left cross axis (positive defined to be right-handed roll about the forward axis). We computed heading as the angle between the +x axis of the world frame (to the right in the video) and the vector formed by projection of the forward cross axis onto the floor plane (counter-clockwise positive).

We levelled the arena before experimentation, and thus the gravity vector was aligned with the world–z axis. At each time point, the dot product of this vector and each of the accelerometer axes was computed, and subtracted from the accelerometer data to leave dynamic acceleration only. The dynamic acceleration waveforms were low-pass filtered at 300 Hz (third-order Butterworth) on all axes.

A phase-based analysis was used to delineate individual steps and to compute average step waveforms and statistics. Phase-based methods are commonly used in many areas of physics and engineering where analyses of periodic time series data are required (Guckenheimer and Holmes, 1997; Kantz and Schreiber, 1999), and are finding increased utility in biology (Revzen and Guckenheimer, 2008; Revzen et al., 2009). These methods analyse periodic time series data to produce an estimate of the instantaneous phase of a periodic signal at each time point. For example, referring to Fig. 2, the highly periodic fore–aft and vertical acceleration traces contribute to computation of a phase value for each time point, a number that rises smoothly from zero to 2π within each step, at a rate corresponding to the local instantaneous rate of oscillation of the waveform. Because each time point in the data has an associated phase value within a single step period, the step can then be divided into small phase bins within which mean values and statistics for each time series can be computed. The phase variable is analogous to normalizing time to a percentage of step, based on identification of features in the data, but because it is computed based on a mixture of time series information, it is more robust than measures based on extremal values, such as peak finding.

We computed the phase variable from raw kinematic data. Fore–aft acceleration, vertical acceleration, the z-component of the forward direction vector and the x-component of the vertical direction vector were detrended through subtraction of their own heavily low-pass filtered waveforms (third-order Butterworth, 16 Hz), and then normalized through mean subtraction and division by their standard deviation. We then summed these waveforms, with the exception of the forward direction vector z-component, which was subtracted, to form a composite time series. This series was Hilbert transformed to produce the phase variable. This linear combination of time series was chosen, through principle components analysis, because it encompassed a significant amount of variation of the data in every trial (44±5.4%, mean ± s.d., N=88 trials). This final phase variable was low-pass filtered using a first-order Butterworth filter at 50 Hz, and points at which the phase variable cycled from 2π to zero were used to count individual steps.

Numerical simulations were carried out using custom-written scripts in Matlab, by integrating equations of motion with the function ode45. Matlab was also used to evaluate the analytical expression predicting the COM acceleration as a function of duty factor given by the model of Alexander (Alexander, 1980).

We analysed a total of 36 trials, six from each of six individuals. With four steps per substrate type per animal per trial, this resulted in a total of 288 analysed steps. An example of an individual accepted trial is given in Fig. 3. The animal traversed the rigid Plexiglas surface, encountered and crossed the compliant membrane, and then regained the rigid surface before ending the trial.

Statistics of step variables are presented in Fig. 4. Cockroaches ran with mean (±s.e.m.) forward speed 37.2±0.6 cm s⁻¹ on the rigid substrate, and mean forward speed 38.0±0.7 cm s⁻¹ on the compliant substrate, values which were not significantly different (repeated-measures ANOVA, P=0.45). The main effects of step number and trial on speed were insignificant, as were all interactions except for that between step number and substrate (P=0.002). Cockroaches exhibited small, opposite trends in speed on the two surfaces, increasing speed with step number on the rigid substrate (3.8% or 1.4 cm s⁻¹) and decreasing speed (1.0% or 0.4 cm s⁻¹) on the compliant substrate. Step frequencies of 24.5±0.6 steps s⁻¹ (mean ± s.e.m.) on the rigid surface versus 24.7±0.4 steps s⁻¹ on the compliant surface were not significantly different (repeated-measures ANOVA, P=0.54), with all other main effects and interactions insignificant. Examination of COM dynamics made it apparent that noticeable changes occurred in the dorsoventral acceleration. We therefore computed and extracted the peak-to-peak dorsoventral acceleration for analysis. The mean (±s.e.m.) peak-to-peak dorsoventral acceleration was 11.5±0.33 m s⁻² on the rigid surface versus 7.7±0.14 m s⁻² on the compliant surface, values which were significantly different (repeated-measures ANOVA, P=0.04). All further main effects and interactions were insignificant.

To understand how the animal's dynamics changed on the compliant substrate, we computed the average time series for several kinematic variables as a function of substrate. Average step time series (means ± s.e.m.) are plotted in Fig. 5. Dorsoventral COM acceleration waveforms exhibited a large single peak with each step, starting at approximately −4 m s⁻² (this is the apex of the COM displacement trajectory, where double support occurred, and the closest the animal came to free fall, which would be −9.81 m s⁻²), and increased to ~+2 to +3 m s⁻² at mid-stance (at the nadir of the COM trajectory, where a single support tripod was generating maximum upward force). On the compliant surface, dorsoventral acceleration was larger at the beginning and end of the step, and smaller through the middle of the step. Rostrocaudal COM acceleration exhibited sinusoidal behaviour, the COM decelerating during the first phase of stance and accelerating in the last. The rostrocaudal COM accelerations were not significantly different between surfaces. These sagittal plane accelerations (a single large peak in the vertical, and sinusoidal, negative then positive phase dynamics in the fore–aft direction) qualitatively resembled those predicted by a SLIP model of a running animal. Lateral COM accelerations for right and left tripod steps displayed rapid transients to acceleration of about 1 m s⁻² during the first 50% of stance, towards the contralateral tripod followed by increasing acceleration towards the ipsilateral tripod. Lateral accelerations did not significantly differ between surfaces.

Forward speed slowed during the initial phase of stance, and then recovered in the second half of stance. We saw larger variation in speed on the elastic surface. The pitch time series showed a minimum just before mid-stance; this corresponded to the animal becoming more vertical, and then returning to a more pitched-back attitude. The animal exhibited the same pitch change on the two surfaces, and again showed greater variation on the elastic surface. The roll time series qualitatively resembled that of lateral acceleration, having symmetric excursions of 25 deg. to the contralateral side of the current tripod (rolling away from the current tripod), during the first 50% of stance, before rolling back to the ipsilateral side. Roll time series were similar on the two surfaces.

We performed numerical simulations of a SLIP model with cockroach-like parameters in order to interpret the observed changes in COM dynamics (Fig. 6). Parameters mass m=2.5 g, elastic spring leg stiffness k=15 N m⁻¹, rest length L=0.024 m, initial horizontal velocity V_x0=0.36 m s⁻¹, initial vertical velocity V_y0=0, landing angle α_SLIP=21.71 deg. and g=9.81 m s⁻² produced symmetric stance phases with sagittal plane accelerations and a step duration similar to that observed in the cockroach running on our rigid surface. We then considered the simplest possible addition of an elastic surface to the model, by placing a spring element of stiffness 10 N m⁻¹, representing our elastic membrane, in series with the leg spring. These compliances in series combine as product over sum to produce an effective leg stiffness of 6 N m⁻¹. Fig. 6B,E (orange dashed lines) illustrates the trajectory of the SLIP with the reduced leg stiffness and initial conditions identical to the rigid substrate. It can be seen that the COM trajectory is no longer symmetric, and that the model falls. It exhibited a lower peak amplitude dorsoventral acceleration, and a step duration much greater than that on the rigid surface.

We found that a symmetric gait can be recovered (Fig. 6C,E, green lines) by increasing the landing angle of the SLIP model (the angle at which the leg touches down, as measured from the surface normal). The dorsoventral acceleration for this set of model parameters reached a greater peak than the rigid surface model, however, and had a stance duration approaching double that of the original, rigid surface model.

Cockroaches run with a short period of double support, such that during the transition between support by one tripod of legs and the other, all six legs briefly contact the surface (schematically shown in Fig. 6D). This feature is not modelled by the SLIP. We therefore used the phenomenological model of Alexander (Alexander, 1980) to predict the effect of changes in the double support period on the vertical acceleration of the COM. The ratio of each tripod's stance duration to the full stride duration is defined as the duty factor, and for cockroaches running at their preferred speed it is slightly greater than 0.5 (Full and Tu, 1990). The model sums sinusoids representing the ground reaction force of each support tripod to predict the overall vertical acceleration of the COM, and varies the period of each sinusoid to simulate longer or shorter double support phases. As duty factor increases, there is a longer double support phase, and the vertical acceleration goes less negative in transition between support tripods (the COM is further from free fall, which would be −1 g) and thus for the net acceleration over a stride to still be zero, the peak acceleration must be lower (Fig. 6D,F). This is achieved in the model by normalizing the amplitude of the summed sinusoids with a multiplicative factor, such that the net vertical acceleration of the COM over a stride remains at zero. With input parameters taken from the experimental data, step duration 44.25 ms and peak-to-peak amplitude of the order of 10 m s⁻², a duty factor of 0.54 accurately reproduced the rigid surface dorsoventral accelerations we have recorded (cf. Fig. 6F and Fig. 5A, blue lines). Increasing the duty factor to 0.57 while holding the other parameters fixed resulted in the red line, now a qualitatively accurate depiction of the accelerations seen on the compliant membrane (cf. Fig. 6F and Fig. 5A, red lines).

Cockroaches, when confronted with an elastic surface having stiffness approximately 2/3 of their leg stiffness, were able to continue forward locomotion at or above their speed on a rigid surface (Fig. 4). Step frequency was unchanged, the body of the animal was less pitched head-up on the elastic surface, and the amplitude of oscillation of the COM acceleration in the dorsoventral axis was smaller on the elastic surface (Fig. 5).

To interpret these changes, we initially turned to the aforementioned simple mechanical model, the SLIP, and asked whether simple changes in one or more model parameters reproduced the changes we observed in our data. We began with this ‘template’ (Full and Koditschek, 1999) because it is the simplest possible model of a bouncing gait; if we found it did not explain our results satisfactorily, we planned to move to a more representative, or ‘anchored’, model. The running cockroach has previously been shown to exhibit SLIP-like COM dynamics during normal running on a rigid surface (Blickhan and Full, 1993; Full and Tu, 1990), a finding which we confirmed upon examination of the shape of our average step dorsoventral and rostrocaudal waveforms (Fig. 5). In the vertical or dorsoventral axis, this comes in the form of a single positive sinusoidal acceleration hump, resembling closely the function sin(t) from t=0→π. In the fore–aft or rostrocaudal axis, we found the expected sinusoidal behaviour with the acceleration resembling –sin(t) for t=0→2π.

The simulation results in Fig. 6 show that a SLIP model tuned to reproduce our rigid surface experimental data does not predict the behaviour of the animal when it is confronted with a compliant surface, and in fact the model becomes unstable. While a stable gait can be recovered through increased landing angle, this compensation resulted in a larger peak dorsoventral COM acceleration and longer step duration. Neither of these changes was seen in our experimental data, wherein the animal exhibited similar step duration and lower peak dorsoventral COM acceleration. Thus, an albeit simple change in SLIP model parameters was not adequate to explain our experimental findings, as was the case in previous studies (McMahon and Cheng, 1990; Ferris et al., 1998). The reduction in COM acceleration (and hence displacement) was also in contrast to the results of the study by Moritz and Farley, who found that humans hopping on very soft elastic surfaces can drive the surface to the extent that their COM displacement is unchanged; although at both hopping frequencies studied a statistically insignificant trend to lower COM displacement was observed (Moritz and Farley, 2005). Similarly, Kerdok and colleagues found no significant change in COM displacement for human runners on surfaces of varied stiffness (Kerdok et al., 2002). For these human runners and hoppers, however, each stance phase was separated by a flight phase, during which time the swing leg(s) is brought forward for the oncoming stance phase. The cockroach studied here did not have a flight phase in between stance phases, and in transition between stance tripods had a brief period of double support.

We then considered how this double-support phase may explain the cockroach's lower COM acceleration using the model of Alexander (Alexander, 1980). With parameters matching our experimental data, this model demonstrated that a slight increase in duty factor, with the required renormalization of the integral of acceleration over time to zero, was enough to predict the changes we observed experimentally (cf. Fig. 6F and Fig. 5A). For the animal that this model described to produce the changes we have computed, it would need to both (1) increase the fraction of time each leg is on the ground during a stride such that the required change in duty factor was met, and (2) develop smaller peak vertical acceleration at mid-stance with each leg (i.e. it would produce the red vertical acceleration curve seen in Fig. 6F, as opposed to the blue curve).

A model actively modulating its duty factor to compensate for the compliant surface may be unnecessary due to the nature of the animal–surface interaction. On inspection of the animal supported by the membrane, it is clear that the membrane deforms to some degree locally about each foot, and that the swing tripod moves forward to contact undeflected membrane. This means that while the cockroach was on the membrane, the COM was effectively closer to the surface. A 2.5 g animal on a 15 N m⁻¹ membrane can be expected to sink 1.6 mm. This is a significant fraction of the ~1 cm hip height of the animal.

Considering the motion of the swing legs, it is clear that for an animal that has ‘sunk into’ the surface, swing legs will hit the surface before they would normally do so on a rigid surface. This will, in turn, result in an increase in duty factor, and this increase would be a consequence of the nature of the mechanical interaction between an organism that is utilizing reciprocating legs and a surface that deforms locally, allowing the COM to sink and protracting legs to make contact with the surface earlier in their swing phase. We hypothesize that the increase in duty factor due to sinking into the surface can ‘automatically’ compensate for the slower loading and force production of the decreased leg–surface system stiffness, simplifying the task of the neural control architecture. While our model showed that on the compliant surface the animal must develop smaller peak vertical acceleration at mid-stance, this could happen as a consequence of different swing leg posture on touchdown.

This potentially self-stabilizing interaction may have important implications for locomotor performance, mechanical stability and neural control. For performance, the work of McMahon and Cheng demonstrated that this mechanism can produce a beneficial increase in stride length (McMahon and Greene, 1979). For mechanical stability, it is important to ask whether a bipedal SLIP model running on a compliant surface has a larger stability basin than the same on a rigid surface (Geyer et al., 2006). For neural control, we turn our attention to whether the nervous system must change its motor commands whilst on the compliant surface, or whether the mechanical interaction ensures that identical instructions will result in successful locomotion. This mechanism has been described in application to the aforementioned study of running humans confronted with a rigid upward step perturbation (Grimmer et al., 2008). The fact that the same leg–surface interaction mechanism may operate to simplify the control task in both bipedal runners with single support confronted with a rigid, increased height obstacle, and hexapedal runners with double support sinking into a soft surface, hints that it could be an important general principle of legged locomotion.

To test mechanistic hypotheses that address the aforementioned questions of control and stability, a model that includes swing leg motion and a hypothetical controller is required. Here we identify such a model through application of a recently formulated model of legged locomotion, the clock-torqued spring-loaded inverted pendulum (CT-SLIP) model, to our results (Seipel and Holmes, 2007). Using this model, we suggest that active control may not be required to produce perturbation responses similar to the cockroach in this simply controlled spring–mass model (Figs 7 and 8). This monopod model has muscle-like torque actuation at its hip and a passive springy leg, which approximates the dynamics of the cockroach, and enables us to determine whether it is possible for feedforward control without active feedback to generate the perturbation recovery. The CT-SLIP can explain the cockroach's perturbation response without any active sensing or feedback. We note, however, that this robust stability could also be achieved by a proprioceptive clock-controlled approach which is similar in dynamics. Further, other methods of state feedback, and finite event feedback at touchdown and lift-off events could also enable robust stability of gait.

Fig. 8B shows the COM path of the CT-SLIP over 15 strides with the membrane interaction approximately from stride five to stride 10. Fig. 8A represents the vertical acceleration normalized by gravity. It clearly shows that the average peak-to-peak acceleration decreases over the membrane, in agreement with our experimental results. The change in ground stiffness was modelled to first-order approximation by lumping the effect of a ‘ground’ spring, in series with the leg spring, into an equivalent springy leg. The purpose of this simple modelling choice was to represent minimally the type of perturbation introduced by the membrane, and to determine whether the model explains the type of response observed experimentally. The response of the model does not change significantly (and does not change in stability type) with small changes in damping or stiffness. Instantaneous spikes in acceleration can occur at touchdown events due to the discontinuities in the hybrid dynamical model, but would not necessarily represent more realistic collision physics. In fact, the inclusion of foot–ground interactions would likely smooth out the forces at transition events.

The antagonist pairing of these feedforward muscles effectively produces a net torque term that can track the error between the leg angle and its equilibrium point dictated by the feedforward signals. However, as in the robotic hexapod RHex (Koditschek et al., 2004), this torque can be generated with a single motor at the hip. But since this motor does not have an antagonist, it must somehow track the leg position and use local proprioceptive control to generate the corrective torque. In RHex this is generated by an encoder which measures the angle of the motor shaft, and a digital proportional-derivative controller which processes the encoder signal and sends an appropriate voltage to the motor. The CT-SLIP model, which closely models the mechanisms of RHex, captures the robust stability exhibited by the robot with only proportional feedback, and in fact with no feedback at all during stance. We hypothesize that other variations of the model, representing different kinds of feedback, can generate robust stability as well. However, we have found that versions of SLIP which have limited mechanical feedback and limited feedforward actuation are not robust, and would not succeed in traversing a membrane. We propose to measure the muscle activation in key extensor and flexor pairs of Blaberus discoidalis in future membrane perturbation experiments to determine whether active changes in neural activation are occurring.

The CT-SLIP model provides a relatively simple hypothetical mechanism as to how the cockroach can successfully traverse the elastic surface while exhibiting the changes in COM acceleration we have measured. Physically, the sudden reduction in leg stiffness causes the animal to compress it further, sinking into the compliant effective spring. The model does not fall, however, because the lowered COM height causes the next swing leg to contact the surface sooner than it would have on a rigid surface. Thus, the more heavily compressed lumped leg spring and lower COM result in increased duty factor, or duration of double support. These changes arise without the need to alter swing leg kinematics, and in fact are the result of maintaining normal swing leg kinematics about a lowered COM.

A consequence of this increased vertical force during double support is that, if we assume that the shape of the vertical force remains constant, then the peak vertical force at mid-stance must be reduced for the average vertical force to still equal body weight. This change does not necessarily require active modulation of stance leg kinematics. We hypothesize the simplest mechanistic explanation for the reduction in mid-stance vertical force as follows. If we assume that the animal has sunk into the more compliant effective leg spring, and consider the force–length relationship for a linear elastic spring F=–k(x(t)–x₀), we have a situation in which k has been reduced, and x(t), composed of a constant offset and an oscillation about that offset, has shifted to a new offset such that F continues to oscillate about body weight. If we simplistically consider that maintaining the same leg kinematics is equivalent to maintaining the same magnitude of oscillation of x(t) about its larger constant offset, then the effect of the reduced k will be to produce a smaller peak F, and potentially the difference we see in peak vertical force at mid-stance. This suggests that it may be possible to maintain leg kinematics on the more compliant surface, both during swing and during stance, yet maintain forward locomotion, and to produce the effects that we have observed in the animal. Furthermore, smaller oscillations in vertical acceleration and hence force may provide a mechanistic explanation for the past findings that insects (Full and Chang, 1995) and humans (Kerdok et al., 2002) consume less metabolic energy when running on a compliant surface, and the performance of human athletes with spring-like prosthetics (Buckley, 2000). In essence, the springy surface does some of the work of redirecting the organism and bouncing it along, relieving to some degree the animals' legs of this responsibility.

The accelerometer backpack we present opens the door to exciting new experiments. Large quantities of continuous COM acceleration data, arriving without modification in the animal's body coordinate system, can be collected. This also allows novel statistical analyses, e.g. quantification of dynamic stability, that require large numbers of contiguous strides (Full et al., 2002; Kang and Dingwell, 2006). Real-time transduction of motion by the accelerometer avails the investigator with perturbations that are phase locked to body dynamics. Through integration of the accelerometer backpack with neural or myographic electrophysiology (Spence et al., 2007; Sponberg et al., 2007), we have a new tool with which to measure the role different parts of the neuromuscular system play in whole body behaviour. Ambitious, systems-level hypotheses of locomotion are arising with the aid of mathematical models such as that of Holmes and colleagues (Holmes et al., 2006). These models extend from ion-channel currents to centre-of-mass dynamics, by way of model neurons, a simple central pattern generator network, and six two-component legs with Hill-type muscle actuation. We hope the accelerometer backpack we present here will yield new, multimodal data with which to confirm or refute these new, integrative neuromechanical hypotheses. The current limitations of the standalone backpack lie in the need to subtract gravitational acceleration to obtain dynamic acceleration, and in signal drift that precludes integration to obtain velocity or displacement over long time scales. Through fusion with other sensors and the continuous refinement of these devices, these restrictions are likely to be overcome in the near future.

The authors would like to thank Dr Simon Sponberg, Tom Libby, Pauline Jennings, Dr Eileen Hebets, the Department of ESPM at UCB, members of the Full laboratory, and the Berkeley biomechanics group for their comments and fruitful discussions. This work was supported by National Geospatial Intelligence Agency (HM1582-06-1-2019) and NSF (EF 0425878) Frontiers in Integrative Biological Research Grants to R.J.F.

 
• a_dyn

  three-element vector of the dynamic acceleration of the accelerometer in the world coordinate system

 
• a_meas

  three-element vector of the measured acceleration signal output by the accelerometer chip

 
• COM

  centre of mass

 
• CT-SLIP

  clock-torqued spring-loaded inverted pendulum

 
• DLT

  direct linear transformation

 
• EMG

  electromyography

 
• F

  force

 
• g

  three-element vector of the acceleration due to gravity, in world coordinates; equal to [0, 0, −9.81]′.

 
• k

  stiffness of elastic spring leg used in SLIP simulations

 
• K

  leg angle proportional feedback gain of the hip torque in the CT-SLIP model

 
• KF

  Kalman filter

 
• L

  rest length of elastic spring leg used in SLIP simulations

 
• LLS

  lateral leg spring

 
• m

  mass of the point mass used in SLIP simulations

 
• MEMs

  microelectromechanical systems

 
• PCB

  printed circuit board

 
• r

  radius from the centre of the calibration turntable at which the backpack was placed

 
• R(t)

  3×3 rotation matrix that transforms from world to accelerometer coordinates; estimated for each time point by the video tracking system

 
• SLIP

  spring-loaded inverted pendulum

 
• t

  time

 
• UKF

  unscented Kalman filter

 
• ν

  linear velocity of the accelerometer backpack on the calibration turntable

 
• V_x0, V_y0

  initial horizontal and vertical velocities used in SLIP model simulations

 
• α_SLIP

  landing angle of elastic spring leg used in SLIP simulations; 0=vertical

 
• α, β, γ

  first, second and third Euler angles, defined in the standard order

 
• τ

  net torque about the hip of the CT-SLIP model

 
• ϕ

  angle of the leg in the CT-SLIP model

 
• ϕ_r

  reference or desired angle of the leg in the CT-SLIP model, set by the clock

 
• ω

  angular velocity of the turntable used in calibrating the accelerometer backpack