Environmental interaction influences muscle activation strategy during sand-swimming in the sandfish lizard Scincus scincus

Successful locomotion in a diversity of environments (Blickhan and Full, 1993; Dudley, 2000; Alexander, 2003; Lauder and Tytell, 2005) emerges from the interplay between neuromechanical (Nishikawa et al., 2007) systems (the combined action of the nervous and mechanical systems) and the physical environment. The central nervous system (CNS), musculoskeletal system and sensory system constitute a feedback system that accepts inputs from internal and external sensors and generates appropriate muscular action to generate forces that act on the environment. At one extreme, locomotion can be generated with a feed-forward control strategy in which activation patterns are dictated by commands from the CNS. This motor control may be insensitive to perturbations (Kandel et al., 2000). In this situation, mechanical feedback can be used to passively stabilize movement with changing external forces (Brown and Loeb, 2000; Spagna et al., 2007; Sponberg and Full, 2008). At another extreme, motor control can be actively modulated by sensory input to generate specific kinematic patterns that are tuned in response to environmental changes (Biewener and Gillis, 1999; Cowan et al., 2006).

The surrounding physical environment plays a significant role in determining animal kinematics, neuromuscular strategy and the efficacy of the locomotion strategy. Many recent studies (Gillis, 1998a; Gillis, 1998b; Biewener and Gillis, 1999; Biewener and Corning, 2001; Gillis and Blob, 2001; Korta et al., 2007; Horner and Jayne, 2008; Tytell et al., 2010; Gidmark et al., 2011) have investigated how these factors change with either changing substrate mechanics or changes in the rate of movement through the material. For example, some animals move between environments like water and land that present different mechanical responses to limb and body movement (Biewener and Gillis, 1999; Biewener and Corning, 2001; Ellerby et al., 2001; Gillis and Blob, 2001). These different force conditions are correlated with different movement strategies. For example, on hard ground, salamanders generate propulsion with limb–ground interactions (Ashley-Ross, 1994) and use standing waves of axial bending (Frolich and Biewener, 1992); however, in water, salamanders generate propulsive thrust forces using traveling wave axial body undulations while limbs are placed along their sides. Lungfish, which commonly move in muddy conditions in the wild, change kinematics and muscle activation when they encounter fluids of varying viscosity (Horner and Jayne, 2008). During steady-state swimming, lungfish increase lateral bending and use a higher amplitude of muscle activity as fluid viscosity increases. Similarly, the American eel (Gillis, 1998a) uses a different pattern of muscle activation and kinematics while swimming in water at faster speeds. Only posterior musculature is employed during slow swimming while anterior muscle groups are recruited as the eel increases its speed.

Modulating neural control strategies with changing conditions can allow an animal to maintain stability, reduce energy, enhance speed or operate at some other locomotion control target. These control targets can be analyzed in more detail by establishing a ‘template’ for a locomotion pattern. Templates are low-order models (Full and Koditschek, 1999) that capture the basic features of a locomotor behavior – for example, an inverted pendulum model for walking (Cavagna et al., 1977) and a spring-loaded inverted pendulum (SLIP) model for running (Full and Tu, 1991; Full and Koditschek, 1999) – and these models have demonstrated predictive ability for stability (Schmitt et al., 2002) and energy savings (Cavagna et al., 1977; Farley and Ko, 1997). Templates can also be used to advance hypotheses concerning neuromechanical control strategies that enable a targeted behavior (Full and Koditschek, 1999).

The biomechanics of the many animals that live on and within granular media, a material that can display both solid and fluid-like features, are relatively unexplored (Norris and Kavanau, 1966; Arnold, 1995; Irschick and Jayne, 1998; Herrel et al., 2011; Korff and McHenry, 2011) and even less is known about the neuromechanical strategies used by animals that move in granular substrates (Jayne and Daggy, 2000). Granular media are defined as collections of particles that interact through dissipative, repulsive contact forces. The animals encounter external resistive force conditions that depend on the location of the movement (above or below the surface), and also on granular media compaction and depth beneath the surface. Compaction is characterized by volume fraction (ϕ; the ratio of the solid volume of the grains to the total occupied volume) (Nedderman, 1992) and is influenced by wind and animal digging, with ϕ ranging from 0.55 to 0.63 in natural dry sands (Dickinson and Ward, 1994).

In the laboratory, ϕ can be systematically varied via an air-fluidized bed (Li et al., 2009; Maladen et al., 2009; Gravish et al., 2010) (Fig. 1B), and provides a controlled method to test the effects of compaction on locomotion. A less compact granular medium (low ϕ, loosely packed or LP) generates smaller forces on a rod intruder moving at constant speed than does a more compact granular medium (high ϕ, closely packed or CP) (Maladen et al., 2009; Gravish et al., 2010) (Fig. 1C,D). Unlike fluids, the average force acting on an object moving at low speed (typically less than 0.5 m s⁻¹) through granular media increases proportionally with depth (Fig. 1C) and is approximately independent of speed (Fig. 1D) (Albert et al., 1999; Maladen et al., 2009). This is because frictional interactions dominate the forces between grains and intruders.

Recent work (Maladen et al., 2009; Maladen et al., 2011b) used high-speed (up to 1000 frames s⁻¹) X-ray imaging to investigate the kinematics of the sandfish (Scincus scincus), a sand-specialist lizard (Mosauer, 1932; Arnold, 1995; Schleich et al., 1996; Attum et al., 2007) native to the Sahara desert (Fig. 1A). The sandfish is subarenaceous, spending most of its time subsurface and moving while fully buried beneath the sand (Mosauer, 1932; Arnold, 1995; Maladen et al., 2009). The X-ray studies showed that the sandfish uses its limbs in a diagonal gait pattern to generate forward propulsion when moving above the surface and, unlike other lizards, the trunk undergoes little undulation (Maladen et al., 2009). As the sandfish dives into the surrounding sand, it executes a stereotyped swimming pattern making it an excellent animal model for studying how the physics of granular media influences movement strategies. During subsurface movement, the sandfish propagates an anterior-to-posterior traveling and approximately sinusoidal wave of axial body undulations to propel itself through the media while its legs remain by its sides (Maladen et al., 2009). Remarkably, the animal uses a similar spatial form (the ratio of wave amplitude, A, to wavelength, λ) of ~0.2, and wave efficiency (η, the ratio of the average forward speed of the animal to the speed of the kinematic traveling wave) of ~0.5 regardless of media ϕ or swimming speed; for a single period sinusoidal wave, A/λ characterizes the shape of the animal. Models reveal that the sandfish propels itself by swimming through a self-generated localized ‘frictional fluid’ within the media (Maladen et al., 2009).

We hypothesized that the sandfish uses a neural strategy that targets these kinematics and we propose that the observed waveform is a template for optimal sand-swimming. In support of this hypothesis, we note that recent theory, numerical simulation and robot models of the sandfish interacting with granular media have revealed that forward velocity is maximized at A/λ≈0.2 (Maladen et al., 2011b). In this paper, to determine whether the sandfish actively modulates its muscle activation strategy to achieve this control target, we measured the kinematics and the axial neuromuscular strategy during above-surface movement and subsurface swimming using high-speed X-ray and visible light imaging with synchronized electromyogram (EMG) recordings of epaxial musculature. Based on above-ground kinematics we predict that there will be little axial muscle activation during walking on the surface. Subsurface, we predict that the sandfish will develop a traveling wave of muscle activation to drive the traveling kinematic wave. Considering the general features of granular forces described above, during subsurface movement we predict that resistive forces in granular media should dominate the muscle activity. Further, if the animal seeks to maintain spatial form in the face of a changing environment, it will have to modify muscle force and potentially its activation pattern. Therefore, as it moves within a frictional fluid we predict a change in EMG with both depth and ϕ but not with undulation frequency.

Sandfish lizards, Scincus scincus (Linnaeus 1758) (Fig. 1A) were purchased from commercial vendors (LLL Reptile, Vista, CA, USA and Ocean Pro Aquatics, Chino Hills, CA, USA). The nine animals used for the EMG and kinematic study had a snout–vent length (SVL) of 8.9±0.4 cm and a mass of 16.6±1.8 g. We did not use sandfish with regenerated tails. Two sandfish were dissected to facilitate understanding of axial musculature (Fig. 2A) and electrode placement (Fig. 2B). The skeletal structure of a live anesthetized sandfish was observed using micro-computed tomography (micro-CT). All animals were housed individually in large containers (21×43×28 cm) filled with sand to a depth of ~15 cm. Animals were fed mealworms coated in a supplemental calcium powder. The holding room was maintained on a 12 h:12 h light:dark cycle and the average temperature was 28°C during the day and 24°C during the night. Animals were provided with a heat lamp and a water dish to facilitate thermoregulation. All experimental procedures were conducted in accordance with the Georgia Institute of Technology IACUC protocol number (A08012) and Radiation Safety protocol (X-272).

Granular media compaction varies in nature and influences the resistive force experienced by a burrowing animal. A few percent increase in ϕ can substantially increase resistive force on an intruder moving through the substrate (Maladen et al., 2009) (Fig. 1D). Further, when an intruder moves through a LP granular medium, drag forces are approximately constant; however, drag forces in CP states can oscillate during the movement through periodic yielding along transiently stable shear planes (Maladen et al., 2009; Gravish et al., 2010). LP and CP states are distinguished by ϕ of less than or greater than a critical volume fraction (Schofield and Wroth, 1968; Gravish et al., 2010), ϕ=0.605 in 0.3 mm diameter glass particles.

To control the compaction of the substrate, an air-fluidized bed (Maladen et al., 2009) (Fig. 1B) was used to create repeatable initial ϕ states. The bed was filled with spherical glass beads that were similar in size (diameter 0.27±0.04 mm) and density (ρ=2.5 g cm⁻³) to particles found in sand dunes (Dickinson and Ward, 1994). An air flow created a fluidized state and a slow decrease in air flow to zero resulted in a LP state (ϕ=0.584±0.013). CP states (ϕ=0.635±0.013) were created by applying a slight air flow below the onset of fluidization and subsequently vibrating the granular media using a motor with an off-axis mass that was attached to the side of the bed. The difference in compaction resulted in an ~50% increase in the average resistive force between LP and CP preparations (Fig. 1C). Air flow was off during all animal trials and ϕ was determined by measuring the height of the granular media before each sandfish swim.

The ambient air temperature during experiments was 34°C. Heat lamps over the sand were turned on at least 1 h prior to recording the first animal trial, resulting in a minimum internal sand temperature of 30°C. The sand bed was fluidized periodically, allowing mixing of the grains, which produced a more uniform temperature distribution within the media. The temperature of the granular media continued to rise during the course of the experiment to a maximum of 38°C. As temperature was measured at the beginning and end of experiments, we used the time stamp on each trial in order to estimate the temperature of the sand and correlate it with speed and EMG (see Results). We estimate that the majority of the experiments (>75%) were conducted between 34 and 38°C. Previous studies on the burial performance of Uma scoparia (Jayne and Daggy, 2000) within granular media at low (~23°C), moderate (32°C) and high (39°C) temperatures revealed that low temperature increased burial time while middle to high temperatures showed similar performance.

Above-surface visible light video was taken at 250 frames s⁻¹ using a high-speed camera (AOS Technologies AG X-PRI, Baden Daettwil, Switzerland) to characterize general features of above-surface movement in the sandfish such as axial bending and speed. Lizards were initially placed in a holding pen (Fig. 1B) separated from the sand-filled container by a gate (4.5 cm wide). Sandpaper bonded to the holding pen floor increased traction during above-surface movement. Two holding pens were used depending on the goal of the study. One had a rectangular area (6.3×18×11 cm) that did not extend over the sand (Maladen et al., 2009) and thus did not interfere with subsurface video recording. In the other holding pen (4.2×11×6.5 cm), the gate was farther from the sand area to induce more above-surface strides before burial took place. In addition, the bottom of this holding pen extended over the sand bed area so that above-surface movement was captured on both the visible light video and subsurface X-ray images. The different holding pens did not affect the subsurface performance of the sandfish. The visible light camera was suspended above the holding pens so that the dorsal view of the sandfish was visible and the midline position could be tracked.

Visualization of subsurface movement in optically opaque material has been accomplished with X-ray (Gaymer, 1971; Gasc et al., 1986; Summers and O'Reilly, 1997) and nuclear magnetic resonance imaging (NMRI) (Baumgartner et al., 2008). In our experiments, sandfish kinematics were recorded using X-ray systems (OEC 9000, Radiological Imaging Systems, Hamburg, PA, USA – top-view recording; and Philips, Andover, MA, USA – side-view recording) coupled to high-speed cameras (Fastcam 1024 PCI, Photron, San Diego, CA, USA and Phantom v210, Vision Research Inc., Wayne, NJ, USA, respectively) which recorded at 250 frames s⁻¹. The top-view X-ray system was set to an energy level of 85 kV at 20 mA and the side view was set to 85 kV at 75 mA. During side-view recordings, the sandfish burial began near the outside edge of the X-ray image intensifier, which had significant distortion. To account for this distortion, a calibration grid with known dimensions was used to transform distorted images using a custom-written image processing code (Matlab, Mathworks, Natick, MA, USA). Because the motion for the top-view recording occurred near the middle of the X-ray image intensifier and camera lens, only a small amount of camera distortion occurred. However, objects closer to the image intensifier appear smaller than objects farther from the image intensifier (i.e. closer to the source). Therefore, in top-view recordings as the animal moved deeper into the material, the sandfish appeared smaller. To correct for this, we calculated a calibration value for depth (where a calibration value of 1 was used for an object at the surface of the sand) by moving an object with known dimensions closer to the image intensifier. We approximated the depth of the sandfish during a particular undulation and used the calibration value to correct calculated distances (i.e. distances and speeds at undulation 2, 3, 4 and 5 were increased by 2%, 5%, 7% and 10%, respectively).

For enhanced contrast, in both visible light and X-ray videos, a maximum of 12 small opaque markers (~1 mm² each, total mass <0.5 g) made from lead were bonded to the midline of the sandfish using cyanoacrylate glue between 0.2 and 1.1 SVL, where the vent is at ~0.75 of the total body length (snout–tail length, STL) of the animal (Fig. 2B, black squares). A maximum of two markers were placed on each limb (one distal and one proximal). These opaque markers were used for tracking during both surface and subsurface movement.

X-ray imaging of the sandfish in the horizontal plane (top-view recording) was used to characterize subsurface kinematic parameters such as spatial form (A/λ), wave efficiency (η), segmental speed (ν_seg) and the phase of the cycle (N=5 animals; Table 1). X-ray imaging from the side view was used to characterize the angle of descent (θ_d) under different packing conditions (N=3 animals; supplementary material Table S1). During subsurface movement, sandfish swam in a 10 cm deep granular bed of glass beads held in a container (with area 22.9×40.6 cm for top-view recordings, and 11.3×40.6 cm for side-view recordings). The bed width was reduced during side-view recordings because X-ray recordings through glass particles beyond 12 cm gave low contrast images. EMG experiments were conducted in the larger bed with the top-view X-ray imaging (Fig. 1B).

To examine neural control commands and changes in muscle actuation we used EMG, which records voltage gradients in muscle. Although it is not directly proportional to tension in muscle, it can be used as a proxy for muscle force development (Loeb and Gans, 1986; Daley and Biewener, 2003). Also, EMG minimally interferes with movement and can be obtained in muscles that cannot be instrumented with tendon buckles (Loeb and Gans, 1986; Biewener et al., 1988).

The sandfish was sedated by gas induction with 5% isoflurane in medical grade oxygen and then anesthetized with an intramuscular injection of ketamine hydrochloride (200 mg kg⁻¹ body mass). Anesthesia was maintained during electrode implantation by continuous administration of 2% isoflurane gas solution. Sandfish were implanted with bipolar hook electrodes (Loeb and Gans, 1986) constructed from two Teflon-insulated, single-stranded, stainless steel fine wires (50 μm diameter, A-M Systems, Carlsborg, WA, USA) which were ~1 m in length. Insulation was removed from the last 0.5 mm of the wire, and the ends of the two wires were separated by ~1 mm. The electrodes were inserted percutaneously using a 26 gauge hypodermic needle into the epaxial musculature located ~4 mm away from the midline (or halfway between the midline and lateral edge, Fig. 2B, red markers) on one side of the sandfish body; the side implanted was chosen randomly between experiments and each side showed similar patterns. For the five animals that were euthanized or did not survive surgery, the implantation sites were verified to be located in the iliocostalis muscle group with longitudinal position varying approximately ±0.05 SVL. However, because animals were not always euthanized following experimentation, the muscle group implanted and longitudinal position could not be verified for all experiments. Five bipolar electrodes were implanted at 0.3, 0.5, 0.7, 0.9 and 1.1 SVL. However, because of electrodes pulling out or shifting, we usually recorded from three to four of these implanted electrode sites that showed a clear EMG burst during burial. Hook ground electrodes were constructed from a single wire and insulation was removed at one end. Five ground electrodes were implanted in the posterior portion of the tail after the 1.1 SVL electrode site (Fig. 2B, blue marker), where adipose tissue was abundant. These electrodes were wound around each bipolar electrode and connected to ground to minimize noise. All wires were strain relieved at multiple dorsal locations by application of cyanoacrylate glue, which adheres strongly to their keratinous skin and minimally affects body movement; 1–3 weeks after implantation, the glue detached from the skin as a result of shedding, and electrodes were easily removed. See Table 1 for a summary of EMG experiments and location of electrodes.

EMG recordings were amplified by 1000 and bandpass filtered between 300 and 1000 Hz using a differential AC amplifier (A-M Systems, Model 1700, Sequim, WA, USA). The signals were recorded at a sampling rate of 2 kHz (National Instruments, USB-6008, Austin, TX, USA) and analyzed through custom-written programs (LabVIEW, National Instruments, Austin, TX, USA). The bandpass filter frequencies and sampling rate were initially determined to reduce movement artifact and increase the number of channels that could be recorded. Analysis with a wider bandpass (10–5000 Hz) and higher sampling rate (10 kHz) revealed power below 300 Hz. Applying a 300–1000 Hz filter and downsampling this EMG signal in software showed EMG intensity decreased by 59±2% but the functional form of the signal and EMG trends did not change, nor were our conclusions affected.

Kinematics were obtained by tracking the opaque markers along the midline of the sandfish body in both above-surface and below-surface video recordings using custom-written software (Matlab) to obtain positions with time. Above-surface sequences were included in the study for analysis only when animals moved farther than two stride lengths forward and used a diagonal gait. Subsurface trials selected for analysis were those in which the mean sinusoidal positions followed a straight path and occurred at least 2 cm away from the edge of the container. The trials in which the sandfish completed less than two undulation cycles or when the sandfish stopped and started swimming while subsurface were not considered. A single period sinusoid was fitted to the sandfish midline markers and those with fits with r²<0.8 were discarded. We accepted undulations in which sandfish completed the full undulation and continued swimming for at least a quarter of a cycle. The sandfish was considered to be subsurface when the animal body and legs were covered by the granular medium. We included EMG bursts in our subsurface analysis with the tail visible (undulation 2) if limbs were no longer in the sprawled posture and appeared to be passive and placed along the sides of the animal. If the legs were used or limbs and body were visible above the surface, the animal was considered to be in the burial phase and these movements were excluded from the subsurface analysis.

The instantaneous forward, lateral and segmental speeds of the animal were calculated between consecutive video frames as the displacement of the 0.5 SVL marker in the forward (ν_for), lateral (ν_lat) and overall direction of motion (ν_seg), respectively, divided by the elapsed time between frames. We used ν_seg in the subsurface analysis to characterize the overall speed and we compared this parameter with EMG characteristics. If the 0.5 SVL marker detached from the animal, the 0.3 SVL marker was used instead and gave similar results. In order to minimize the effects from the contralateral muscle activation, the mean segmental speed () was calculated as the average segmental speed for approximately half of the total cycle (during the time the muscle at the 0.5 SVL location was shortening). Assuming the sandfish is moving in a perfect sine wave, when ν_lat>0 there is a change in the configuration of that segment from maximally convex to maximally concave; thus, muscle shortening coincides with positive lateral speed.

As in Maladen et al. (Maladen et al., 2009), the midline positions fit well to a single period sinusoidal wave y=Asin[2π/λ(x+ν_wt)] during subsurface swimming, where A is the amplitude, λ is the wavelength, ν_w is the speed of the traveling kinematic wave, x is the position in the forward direction and t is time. A/λ and wave efficiency, η=ν_x/ν_w, were calculated for each undulation cycle. ν_x is the average forward speed and ν_w=fλ, where f is the undulation frequency. For top-view X-ray trials, where a direct measurement of depth by side-view X-ray was lacking, the estimated depth was binned according to undulation number using both above-surface video with top-view X-ray images. Undulation 1 was defined as when the head of the animal was covered by granular media but hindlimbs remained visible. Undulation 2 was defined as when only part of the tail was visible and the rest of the body and legs were subsurface. Undulations 3, 4 and 5 were defined as the first, second and third burst when the sandfish was undulating fully subsurface. Only EMG bursts during undulations 2–5 were considered in the depth and ϕ analysis, while undulation 1 was ignored because a significant proportion of the body was above the surface and use of hindlimbs was evident. In the side-view X-ray recordings, depth was measured directly as the vertical distance from the surface of the media to the 0.5 SVL marker position. θ_d was measured for side-view recordings as the angle between horizontal and the best linear fit to the positions of all tracked dorsal markers during the subsurface descent.

EMG signals were analyzed using SpinalMOD, a custom-written Matlab program (Gozal, 2010). For burst detection, the EMG signal was rectified and filtered using a second-order Chebyshev filter (Fig. 2C, red trace) with a stopband ripple of 10 dB and stopband-edge frequency of 30 Hz. The ‘filtfilt’ function in Matlab was used to remove phase distortion. A threshold equal to the mean of this filtered-amplified EMG trace was set to detect EMG burst. Burst onset was defined as the time when the filtered EMG signal exceeded the threshold and afterwards remained above it for a minimum of 0.04 s (Fig. 2C, magenta marker). EMG burst offset was defined as when the filtered EMG signal became lower than the threshold and remained below for at least 0.08 s (Gozal, 2010) (Hochman et al., 2012) (Fig. 2C, green marker). Threshold and timing parameters were changed when necessary by visual observation because of the occasional occurrence of spikes that did not constitute an EMG burst or missed detection of EMG onset. EMG duration was calculated as the time between the onset and offset of the EMG burst (Fig. 2C), and the EMG duty cycle was determined as the ratio of EMG duration to EMG period (time between sequential EMG bursts). The frequency of bursting was defined as the inverse of the EMG period.

The rectified integrated area was found for each EMG burst using the original EMG signal (filtered once in hardware). Because the area was influenced by burst duration and swimming speed, the EMG intensity (I), defined as the ratio of the EMG integrated area to the EMG duration, was calculated to characterize the differences in activation between different ϕ, undulation numbers and speeds. To reduce variation associated with averaging intensities across multiple recording sites, EMG intensities were only considered for the 0.5 SVL location. Because of the differences in electrode construction, I was compared across individuals only after standardization (see Gillis, 1998a). After observing correlation between undulation number and I, a standardization was performed by using the best-fit line for I versus undulation number for a single implantation site. We used the calculated I at undulation 3 to normalize all other EMG intensities at that site. This undulation number was used because all animals tested consistently undulated at least 3 times during burial. We refer to the normalized EMG signal as the relative EMG intensity, I_rel.

The activation timing was compared with the estimated muscle strain cycle at a segment by calculating the angle between adjacent kinematic markers on the dorsal midline. Using this sinusoidal trajectory of angle with time, the normalized muscle length (L_M) and total cycle duration were estimated. L_M of 0 represents when the muscle is at resting length and the local segments (i.e. axial body positions at 0.1 SVL increments) are aligned. L_M of 1 is when the segment is at maximal convexity and L_M of −1 is maximal concavity. Therefore, moving from L_M=1 to −1 (phase between 90 and 270 deg) represents muscle shortening. The difference between the time at which L_M was maximal and the EMG onset time was divided by the cycle period to calculate relative phasing, which will be referred to as neuromechanical phase lags. Comparing EMG onset with the estimated muscle strain cycle is a common measure among undulatory swimmers (Gillis, 1998a), and it can be used as a means of comparing muscle activation timing in water to activation timing in granular media. We calculated the ratio of the average wave speed of activation to the average wave speed of curvature (ν_EMG/ν_c). ν_EMG was calculated by fitting a line between the longitudinal distance between EMG recording locations along the body and EMG onset time. ν_c was calculated by fitting a line between distance and time at which maximal convexity occurs at those locations. Note that ν_w, the wave of lateral displacement, is 0.68±0.18 of ν_c calculated from curvature speed. A ratio of ν_w/ν_c<1 is expected as the wave of curvature must travel faster than the wave of lateral displacement. This is because the curvature propagates along the entire body (STL) while the lateral wave propagates a distance λ, which is always less than the STL for an inextensible animal. Theoretically, for A/λ of 0.22 this ratio should equal to 0.73.

We used a previously developed granular resistive force theory (RFT) model (see Maladen et al., 2009; Maladen et al., 2011b) to explain some features of sand-swimming. The RFT model used a simplified model of the sandfish kinematics and empirical force laws of granular media as inputs. The modeled sandfish was assumed to have a uniform body and moved in 0.3 mm diameter glass particles. The body of the model sandfish was divided into infinitesimal segments. By prescribing the different orientation of segments as a function of time, waves with different A/λ could be produced. For each segment, forces were calculated as functions of segment velocity and orientation using empirical force laws [see Maladen et al. for functional forms of the force laws (Maladen et al., 2009; Maladen et al., 2011b)]. The forward speed of the body was then obtained from the balance of thrust and drag in the forward direction. The weight-specific mechanical cost of transport (CoT) was calculated by summing the power output of all segments to the external media and dividing by the forward speed and sandfish body weight. The RFT model generated predictions of the dependence of η, mechanical CoT and swimming speed on A/λ for a given undulation frequency, f. By increasing the drag force on the head, the RFT allowed simulation of sandfish movement under different drag conditions, important to bound the effects of drag induced by the EMG wires.

We used analysis of variance (ANOVA) to determine whether spatial form (A/λ) varied with ϕ (LP or CP), undulation number (2, 3, 4) or individual. Next, three multivariate analyses were performed to resolve which factors affected the kinematic variable wave efficiency, η. In the first, analysis of covariance (ANCOVA) was used with A/λ as the continuous variable, ϕ (LP or CP) and undulation number (2, 3, 4) as fixed variables and individual as the random variable. We used ANOVA to test whether η and A/λ in our data set differed from previous data (Maladen et al., 2009) and included ϕ as a factor. To test the effect of wire drag on the response variable, η, we used ANOVA and included the factors of A/λ (continuous), individual and whether the trial did or did not have wires glued to the animal's dorsal surface (fixed).

The temperature effect on animal speed and EMG intensity was characterized by using two ANCOVA analyses. For the first, we tested whether the factors of estimated temperature (continuous), undulation number and individual influenced animal speed. We did this for both LP and CP trials separately. For the second, estimated temperature, speed, ϕ and individual were tested against EMG intensity. In this analysis, we examined statistical results at each undulation number.

A paired t-test was used to determine whether EMG intensity during walking above surface was statistically different from that when remaining still. To see how the parameters of undulation number, ϕ, speed and individual affected EMG intensity during subsurface movement we performed two ANCOVA. In the first we included the interaction term of speed×ϕ and grouped by undulation number. In the second we included the interaction term of speed×undulation number and grouped by ϕ. We repeated the first test but also included the continuous variable A/λ to see whether this affected results. To understand timing characteristics, we calculated the phase lag between EMG onset relative to the time at which that site was maximally convex and the EMG offset relative to maximal concavity (N=4 animals). We used ANCOVA to test whether the response variable, phase (onset and offset), changed with the position along the body (fixed), ϕ and speed. Only cycles during undulation 3 were included in the analysis. Statistical tests were conducted for each of the four individuals separately. Animal 5 (see Table 1) was excluded from the statistical test because of the overall low number of trials collected and, thus, low degrees of freedom. We tested whether the duty cycle of the EMG changed with the factors of undulation number, position along the body and individual using ANOVA. We conducted this test for both LP and CP states. Lastly, ANOVA was used to determine whether undulation number, ϕ or individual affected the ratio between the speed of the EMG wave to the speed of curvature (ν_EMG/ν_c).

In the side-view experiments, we tested whether ϕ and individual influenced the angle of burial using ANOVA. For all above tests, ϕ and undulation number were fixed variables, individual was a random variable and speed was a continuous variable. Statistics were performed using the software package JMP (SAS Institute, Cary, NC, USA). The results were considered significant if the P-value was smaller than 0.05. Values are shown as means ± s.d.

A micro-CT scan of a single sandfish revealed 26 vertebrae in the trunk and 13 anterior caudal vertebrae in the tail; individual posterior caudal vertebrae could not be distinguished. Dissection revealed qualitatively similar muscle morphology to that described for Iguana iguana (Carrier, 1990; Ritter, 1996). The iliocostalis musculature (targeted for implantation) is located on the dorso-lateral portion of the trunk (Fig. 2A).

Fig. 3A shows the qualitative difference in trunk trajectories between above-surface walking (green), burial (pink) and subsurface swimming (blue). On a hard surface before burial, the sandfish moved with a diagonal gait (each front limb moved in synchrony with the contralateral hindlimb) and showed little lateral bending of the body (when speed <3 SVL s⁻¹) (Fig. 3A,D, ‘Above’ region). For a sandfish that was unconstrained (i.e. moved at least 5 cm away from walls and had no wires attached to its body), the maximum lateral deviation of an axial segment from its mean position during walking on a hard surface was 0.61±0.17 cm (n=6). The sandfish with EMG wires showed a slightly greater deviation due to a curved tail from the wire attachment and also tended to turn more as a result of container narrowing. During walking above surface on hard ground, most sandfish moved at slow speeds (relative to subsurface movement) with a median average speed of 0.49 SVL s⁻¹ (n=17 trials, N=4 animals) and a mean speed of 1.25 SVL s⁻¹ (n=17 trials, N=4 animals). All trials had an average speed less than 3.0 SVL s⁻¹ except for one, which had an average speed of 5.8 SVL s⁻¹. For this outlier, lateral bending was larger. After reaching the granular bed, the sandfish paused for 0.8±0.7 s on the surface of the media before beginning burial.

The forward and lateral velocities are shown for a representative trial (Fig. 3F). Lateral velocity was small and the mean forward walking speed was 2.7 cm s⁻¹ (or 0.3 SVL s⁻¹). Forward velocity decreased to zero as the sandfish paused for 0.6 s before entry. The sandfish then tilted its head downwards and thrust itself into the media using its limbs.

During burial, the sandfish also increased the curvature of its trunk (Fig. 3B,D, ‘Burial’ region). It accelerated from rest and buried itself in 1.0±0.3 s (n=17, N=4). In accord with previous findings (Maladen et al., 2009), the hindlimbs were in a sprawled position during above-surface movement and most of burial (distance between the hindlimb distal markers on ankles was 2.8 cm, Fig. 3G). After half of the body was buried, the hindlimbs simultaneously adducted and the distance between limbs decreased. During subsurface locomotion the distance between limbs oscillated around a mean distance of 1.9 cm (Fig. 3G). We attribute the changing limb distance to external forces, as the limbs were closest to the midline when the body location was maximally convex. Forelimbs show a similar trend (Maladen et al., 2009).

During subsurface movement the animal propagated a traveling wave of undulation down its body to move forward (Fig. 3C,D, ‘Subsurface’ region). Sandfish segmental speed was between 1 and 3.5 SVL s⁻¹ during subsurface movement and higher than forward walking speed. Sandfish decelerated during subsurface movement by 0.64±0.4 SVL s⁻¹ (n=21, N=5) for each undulation cycle. Most trials were conducted in a temperature range between 34 and 38°C. However, a few trials were recorded at temperatures as low as 30°C. We found no significant effects of temperature on animal speed.

Average spatial form, A/λ, for animals with EMG electrodes was 0.18±0.04. In accord with our hypothesis, spatial form did not vary with ϕ (ANOVA, F_1,42=2.3, P=0.14), or undulation number (F_3,42=1.4, P=0.26) (supplementary material Fig. S1). However, we did find that spatial form varied between animals tested (F_4,42=3.7, P=0.01). Similarly, wave efficiency, η, was independent of ϕ (ANCOVA, F_1,24=1.4, P=0.25) and undulation number (F_2,24=0.9, P=0.42) (supplementary material Fig. S1), but varied among different animals tested (ANCOVA, F_5,24=5.4, P<0.01). In addition, η was positively correlated with A/λ (ANCOVA, F_1,24=6.5, P=0.02) which is in accord with RFT model predictions (Maladen et al., 2011b) (see supplementary material Fig. S2). All data were in a range where η was sensitive to A/λ (e.g. a 30% increase in A/λ led to a 75% increase in η).

Overall, subsurface sandfish kinematic results were in accord with those found previously (Maladen et al., 2009) and suggested that the EMG wires did not impede the performance of the sandfish. η and A/λ were compared with previously collected data (Maladen et al., 2009) in which no electrodes were inserted into the animal. η and A/λ were higher for the animals in the previous study (η=0.51±0.10, A/λ=0.21±0.04) than for the animals with implanted electrodes (η=0.36±0.08, P<0.01, and A/λ=0.18±0.04, P<0.01) (supplementary material Fig. S2). To test whether drag from the electrode wires affected η, we tested sandfish with (n=11 trials, N=2 animals) and without (n=15 trials, N=4 animals) electrode wires attached to their dorsal surface (no electrodes were implanted). For the animals used during this experiment, we found that η did not significantly vary with individual. We found that there were small but significant differences in η with and without wires (ANCOVA, F_1,20=6.6, P=0.02) but that all η were within the range of η for EMG-implanted sandfish. In addition, best-fit slopes between η and A/λ were similar for sandfish with and without wires (see supplementary material Fig. S1). These experiments were conducted in LP treatments, and undulation number was not characterized because previous findings showed that it did not influence η or A/λ. The increased drag resulted in a lower η (~19% decrease) for a similar A/λ. This decrease in η was in accord with the RFT prediction (supplementary material Fig. S2) that the extra drag from the wires (estimated to be 40% of the head drag by dragging the wires through the granular medium) leads to an ~15% decrease in η in the range of A/λ used by the sandfish (A/λ≈0.2). The animals implanted with EMG wires used in this study displayed a larger range of η than in the previous study; our animals also operated at a lower average η (supplementary material Fig. S2).

EMG recordings were taken from five different sandfish. A total of 33 subsurface top-view trials with EMG were considered (14 CP, 19 LP) from all recorded trials. Within the top-view trials, 54 fully subsurface, straight undulations were recorded and analyzed. We found no significant effects of temperature on EMG intensity.

During above-surface movement, no muscle activity was detected (Fig. 3E, ‘Above’ region). I during above-surface movement (4.1±1.3 μV) did not differ from the noise produced when the wires were moved but the animal remained stationary (5.4±4.8 μV, P>0.05). During the first half of burial, the tracked segments moved forward with small lateral velocity and minimal activation (Fig. 3E,F, first half of ‘Burial’ region). As the animal entered the medium, lateral oscillation increased and EMG bursts were detected (Fig. 3E,F, second half of ‘Burial’ region). During subsurface movement, the sandfish produced an anterior to posterior traveling wave of muscle activation, which coincided with the traveling wave of curvature. The sandfish body shape was similar at the onset of each EMG burst. Muscle activity onset (colored circles) is plotted against lateral velocity at the 0.5 SVL segment in Fig. 3F (‘Burial’ and ‘Subsurface’ region). Near the 0.5 SVL EMG onset (cyan circle), the lateral velocity of the 0.5 SVL segment became positive; a positive lateral velocity indicates muscle shortening, assuming the sandfish is moving in a sine wave.

Typically, EMG bursts became larger as the sandfish descended into the granular medium; this can be seen in a representative trial in Fig. 4A. EMG intensity, I, increased with undulation number for most animals (the slope of I versus undulation number was significantly higher than zero, P<0.05, for the individuals shown; Fig. 4B–D). For combined data (N=5 animals), in which speed, undulation number and individual are factors and the data are grouped by ϕ, we found that relative intensity (I_rel) increased significantly with undulation number (F_2,18=3.7, P=0.04 for LP and F_3,13=5.2 and P=0.01 for CP, Fig. 4E). I_rel increased by 0.26 per undulation in CP and by 0.2 per undulation in LP, but these slopes were not significantly different. Individual, speed and interactions did not significantly change I_rel.

The mean intensity of CP data was higher on average than LP data for every undulation number in the combined data set. These differences were small, with a 13% increase in I_rel from LP to CP during undulation 2 and 3 and a 19% increase during undulation 4. In a test in which speed, individual and ϕ were factors, we found that this intensity was not significantly different between different ϕ at each undulation number (F_1,5=0.2, P=0.64 for undulation 2, F_1,19=0.03, P=0.85 for undulation 3 and F_1,4=5.3, P=0.08 for undulation 4). As in the previous analysis, individual was an insignificant factor influencing I_rel.

Mean speed () also did not influence I_rel (Fig. 5A) for undulation 2 and 4 (F_1,5=1.2, P=0.33 and F_1,4=6.9, P=0.06, respectively) but was calculated as a significant factor for undulation 3 (F_1,19=7.0, P=0.02). However, for undulation 3 the slope of the linear fit between I_rel and speed was small (0.15 s⁻¹) with a low coefficient of determination (r²=0.13). Including spatial form in the analysis did not affect these results. was strongly correlated with undulation number, with a smaller for a larger undulation number (F_3,16=10.8, P<0.01 for CP and F_2,20=15.3, P<0.01 for LP). These trends were observed in each individual (Fig. 5B–D). In accord with previous work (Maladen, et al., 2009), , and increased with undulation frequency.

EMG onset at 0.3 SVL began after the muscle began to shorten (Fig. 6C, magenta bar). The more posterior muscles (0.5–0.9 SVL) activated earlier relative to their muscle strain cycle (i.e. while the muscle was lengthening). The neuromechanical phase lags between EMG onset relative to curvature increased in the posterior direction for all animals statistically tested (N=4; P<0.01). In two animals, there was a significant effect of both speed and ϕ of onset phasing (animal 3: F_1,7=7.8, P=0.03 for speed and F_1,7=10.2, P=0.02 for ϕ; animal 4: F_1,7=7.7, P=0.03 for speed and F_1,7=6.3, P=0.04 for ϕ). Tests show that, on average, phase shifted to the right (EMG onset occurs later in the strain cycle) when the animal moved at higher speeds, while there was a shift to the left when the animal moved through the higher ϕ. However, this was not significant for the other two animals, which had a higher number of overall trials (animal 1: F_1,13=0.2, P=0.69 for speed and F_1,13=1.0, P=0.33 for ϕ; animal 2: F_1,25=4.2, P=0.05 for speed and F_1,25=0.3, P=0.59 for ϕ). In accord with onset timing, EMG offset phase lags relative to curvature statistically increased in the posterior direction for animals 1, 2 and 3 (P<0.05), but not for animal 4. Speed and ϕ were insignificant factors affecting offset phase lags for all animals tested except for animal 3 (F_1,5=31.0, P<0.01 for speed and F_1,5=33.8, P<0.01 for ϕ). As a result of the increasing neuromechanical phase lags along the body, the average wave speed of muscle activation (ν_EMG) was greater than the average speed of the kinematic wave of bending (ν_c) for almost all trials (Fig. 6B), with a median ratio of ν_EMG/ν_c=1.24. We report the median because the distribution is positively skewed. The lower quartile of the ratio is 1.13 and the upper quartile is 1.37.

EMG activity was detected for the majority of muscle shortening in all implantation sites sampled. The EMG duty cycle (EMG duration/EMG period) decreased in the posterior direction (F_3,30=4.2, P=0.01 for CP and F_3,34=4.2 and P=0.03 for LP) from 0.52±0.07 at 0.3 SVL to 0.40±0.06 at 0.9 SVL (Fig. 6D) and showed no dependence on undulation number or individual. Burst duration was inversely correlated with undulation frequency but showed no correlation with I_rel, which supported the finding that I_rel was independent of speed.

Subsurface side-view recordings revealed that markers near the dorsal surface of the sandfish followed a nearly linear trajectory (best linear fits had a mean r²=0.87±0.08, Fig. 7A,B). The sandfish undulated within a plane with angle θ_d with respect to the horizontal. The animal did not pitch around its center of mass during descent into the media. The sandfish dived into the medium at a shallower θ_d (Fig. 7C) in CP media (18.9±4.3 deg, N=4 animals, n=18) than in LP media (26.1±5 deg, N=4 animals, n=18, F_1,28=18.0, P<0.01). The angle also varied significantly with individual (F_2,28=5.9, P<0.01). The average depth between the surface level and 0.5 SVL dorsal marker at which the sandfish first stopped moving was significantly different between LP (3.01±0.93 cm) and CP (1.85±0.54 cm) media (Fig. 7D, P<0.01). Angles and depths reported previously (Maladen et al., 2009) (θ_d=22.2±3.7 deg, depth=2.1±0.5 cm for both LP and CP sand) fall within the ranges observed in this study. In both CP and LP media, the sandfish could dive to the bottom of the sand bed.

Above the surface, the sandfish exhibited minimal body undulation and no epaxial muscle activation; thrust was achieved through leg–sand interactions. The sandfish has prominent limbs (~20% of total body length), suggesting that these appendages are important to sandfish movement strategies. In the sandfish, epaxial muscles are likely involved in lateral bending because epaxial activation coincided with large body undulations during burial and subsurface swimming. The use of epaxial musculature for bending in the sandfish resembles the case in other taxa [e.g. amphibians (Frolich and Biewener, 1992) and snakes (Jayne, 1988)] that use epaxial muscle for lateral bending. There is some evidence that epaxial muscle may be used during high-speed (>3 SVL s⁻¹) above-surface running, but further investigation is needed. More than half of the above-surface trials recorded had an average speed of less than 0.5 SVL s⁻¹ and all trials except one had average speeds <3 SVL s⁻¹. For the one trial that had an average above-surface speed of 5.8 SVL s⁻¹, the animal exhibited some lateral bending and EMG activity was detected beyond the level of noise (I=0.036 mV). However, this intensity was much less than that recorded during subsurface swimming at lower speeds. Below 3 SVL s⁻¹ sandfish did not show significant epaxial muscle activity.

Other studies have shown epaxial muscle activity in lizards [Iguana iguana (Ritter, 1996) and Varanus salvator (Ritter, 1995)] during walking at speeds less than 1.5 SVL s⁻¹. In these studies, the epaxial muscles were found to stabilize the trunk as they were active during the rear foot stance phase during normal walking. The difference in epaxial muscle use could arise from the difference in anatomy and body posture of the sandfish, which may lead to variation in ground reaction forces. The sandfish uses a sprawled gait with its body close to and occasionally touching the ground. The artificial setting of hard ground may differ from the natural above-ground movement that would occur on soft sand. Also, we could not confirm whether the hypaxial muscles were also being employed during movement because this layer of muscle is smaller than the epaxial musculature and difficult to implant with our current EMG electrodes. Future studies that investigate axial muscle use during running on sand over a range of speeds would help elucidate the role of the sandfish axial musculature in stabilization and/or bending.

Jayne and Daggy hypothesized that movement through loose, fluid-like sand could elicit an ancestral traveling wave axial activation pattern in lizards (Jayne and Daggy, 2000). In accord with this hypothesis, we found that the sandfish achieved forward propulsion by propagating an anterior-to-posterior traveling wave of activation along its epaxial musculature during burial and subsurface movement. This wave of muscle activation correlated with the kinematic wave of body undulation. Remarkably, we found that A/λ did not change with undulation number (depth), speed or ϕ. We previously hypothesized that the animal targeted this particular kinematic waveform. We thus predict that to maintain these kinematics as environmental conditions change, the activation strategy must change to accommodate the resistive demands from the environment. We argue that for sandfish-scale animals moving in highly resistive environments like granular media, external resistance dominates while effects of body and material inertia are small; previous studies (Maladen et al., 2009; Ding et al., 2012) estimated that inertial forces at sandfish undulation frequencies were less than 10% of frictional forces. In addition, internal forces from deformation of elastic elements (like tendons) and friction from joint movement will influence muscle activation. However, with similar kinematics we assume that changes in internal forces are small compared with the change in external force due to differing depth and ϕ (Y.D., S.S.S., K. Wiesenfeld and D.I.G., unpublished). Thus, assuming muscle force scales with activation (a point we will return to) the change in activation should correlate with changes in external resistance when kinematics are similar. We therefore predict that resistance forces should give a qualitative indication of the intensity and timing of muscle activation, and we now discuss how depth, speed and volume fraction affect activation strategy.

From the top-view (dorsal) video and EMG measurements, we concluded that as the sandfish swam into the granular medium, I_rel increased with depth (characterized through undulation number) while A/λ did not. Rod drag experiments (Maladen et al., 2009) revealed that for horizontal movement through a granular medium of a cylindrical steel intruder with coefficient of friction comparable to that of the sandfish and at comparable depths, drag force increased approximately linearly with depth by 1.32±0.24 N cm⁻¹ in CP media and by 0.86±0.19 N cm⁻¹ in LP media (Fig. 1C). Also, resistive force in a granular medium and average muscle activation intensity were positively correlated. We hypothesize that for the sandfish, increased I generates larger muscle forces to accommodate the increase in granular force with depth and allows the animal to maintain similar kinematics. The increase in resistive force between estimated depths at undulation 2 and undulation 4 was 75% and 67% for CP and LP, respectively; the increase in I_rel was 63% and 54% for CP and LP, respectively. In contrast, during movement through water, resistance force remains approximately constant for changes in depth, and thus small changes in depth should not change muscle activation strategy for aquatic swimmers. Depth has only been considered for aquatic swimmers during large changes in depth (~800 m) such as in diving sperm whales (Miller et al., 2004).

Our results showing a change in EMG with depth are in accord with other studies of burrowing in soil (Paggett et al., 1998). For example, larval lamprey display an increase in EMG burst amplitude (difference between the maximum EMG burst height and the baseline) when moving from steady-state swimming to burrowing and then during initial stages of burial compared with the final stage. Changes in kinematics also occur between these three states (swimming, initial and final burrowing), so the increase cannot be solely attributed to environmental differences. However, these studies suggest depth of movement is a crucial factor for locomotion in a granular environment. It is interesting to contrast the data for sandfish with a recent study of razor clam digging. Winter and Hosoi recognized that increased depth leads to increased force in wet granular media and suggested that the Atlantic razor clam (Ensis directus) uses a movement strategy that diminishes this force dependence by fluidizing the media surrounding its shell (Winter and Hosoi, 2011).

Although A/λ did not change, we found that the sandfish decreased its undulation frequency as it moved deeper into the material. The animal slowed by 0.64±0.4 SVL s⁻¹ (n=21, N=5) for each undulation cycle. This deceleration is expected to have little effect on resistive force as inertia is negligible. A similar phenomenon was found in the nematode; undulation frequency decreased with increasing external resistance as a result of increased viscosity (Korta et al., 2007; Fang-Yen et al., 2010). This change in frequency may act to limit muscle power expenditure. Korta and colleagues found that mechanosensory input may act to regulate this temporal frequency and allow gait adaptation in response to different surroundings (Korta et al., 2007). Studies of our sandfish models (Maladen et al., 2011a; Maladen et al., 2011b; Maladen et al., 2011c) could be used to determine the extent of muscle power reduction as a result of decreasing frequency with depth for a sand-swimmer.

Like other swimmers, the sandfish increased swimming speed by increasing undulation frequency. However, unlike animals moving in water (Gillis, 1998a), I for the sandfish was insensitive to segmental speed (Fig. 5). We can understand this with the following argument: if the force to move is dominated by body inertia, then we expect that the force should scale like the undulation frequency squared. If, however, the force is dominated by the force to move the environment, then the force should scale like the environmental forces do with speed. Because during sand-swimming inertial forces are small and the granular resistive forces are insensitive to speed, muscle force in a granular medium can be predicted by environmental interaction forces. We note that sandfish decelerate as they move deeper into granular media (Fig. 5); this deceleration is not expected to affect resistive forces from the surrounding media because inertia is small. However, these accelerations may affect the EMG–muscle force relationship.

The above arguments rely on the assumption that muscle activation and force scale proportionally. This might not be true because muscle contractions with differing speed, muscle lengths or activation history may generate different forces for the same activation intensity. When moving at a higher speed, it is possible that the change in sandfish I is due to inherent muscle mechanics (e.g. contraction speed). However, our findings, which show that muscle activation intensities are insensitive to speed, suggest that these effects may be small compared with the variation between trials. In contrast, in a study of the American eel Anguilla rostrata (Gillis, 1998a), which moves with speeds comparable to those of the sandfish, but in a fluid where resistive force increases with speed, muscle activation intensity increased and kinematics changed with increasing speed. This change could result from the changing resistive force imposed on the eel and/or from the muscle mechanics, which are speed dependent [e.g. Hill muscle model (Hill, 1938)].

Surprisingly, although changes of a few percent in ϕ led to relatively large charges in drag force, we did not observe a correspondingly large increase in I_rel for swimming in LP relative to CP. A few characteristics of the animal and the response of the granular medium could explain the smaller than expected changes in activation.

First, the similar activation level could result from the body interacting with disturbed material as the animal proceeds through the granular medium. As the sandfish swims, its anterior body elements shear nearby grains, which causes a CP medium to dilate and a LP medium to consolidate, both evolving toward an intermediate critical packing state (ϕ≈0.605) (Gravish et al., 2010; Umbanhowar and Goldman, 2010). The local volume fraction surrounding the bulk of the swimming animal might be similar in the two preparations despite different ϕ. Consequently, the resistive force could be similar for sandfish swimming in states with different initial ϕ.

Another explanation of the similar activation level could be the different angles of descent between LP (~26 deg) and CP (~19 deg) states. A given undulation number in LP occurs at slightly deeper depths than in CP because of the steeper descent angle in the former; for example, during undulation 3 it is estimated that LP trials are 0.34 cm deeper than CP trials as depth was not directly measured in top-view recordings. The confounding effects of depth may have shifted our results, increasing the average intensity of the LP trials at each undulation number. Thus, there may be a greater I_rel difference between LP and CP than reported because depth was not directly measured. Future work that uses models of both the sandfish and granular media (Rapaport, 2004) could better elucidate how internal forces change with changing parameters such as θ_d, ϕ and depth of sand-swimming. In addition, models could be used to analyze how ϕ changes within the media as a consequence of sandfish movement.

In the swimming sandfish, the average speed of the traveling wave of activation was faster than the average speed of the curvature wave (median of ν_EMG/ν_c≈1.2). For undulatory swimmers moving in true fluids, this difference in average speeds has also been observed and the ratio is 2.4 for trout (Williams et al., 1989) and 1.4 for lamprey (Williams et al., 1989; McMillen et al., 2008). The speed differences occur as a result of changing EMG timing relative to curvature. The phase lag between EMG onset and maximum curvature (neuromechanical phase lag) increases along the length of the animal. Many recent studies have reproduced the neuromechanical phase lags in simulation. For example, Tytell and colleagues simulated an undulatory swimmer and varied undulation frequency, muscle stiffness and fluid viscosity (Tytell et al., 2010). They found that increasing neuromechanical phase lags down the body resulted when average muscle forces were small compared with average fluid forces. A similar phenomenon could explain the phase lags found in sandfish. Other work has suggested that geometry and passive visco-elasticity may affect the extent of the changing phase lags (McMillen et al., 2008; Fang-Yen et al., 2010). For the sandfish, muscle shortening begins in anterior locations before activation is detected. We have found similar results with our RFT model, which prescribes a kinematic waveform similar to the sandfish and couples it to a theoretical model of the granular media. We can use the RFT to predict required torque. We find that bending can start while the torque is produced in the opposite direction in anterior segments because of external forces exerted on the body. A more detailed description will be provided in a future paper which compares EMG timing in sandfish to torque generation in a model undulatory swimmer moving through granular media (Y.D., S.S.S., K. Wiesenfeld and D.I.G., unpublished). Currently, our results are inconclusive regarding the effect of speed and resistive force on activation timing. The RFT model could be used to analyze the effect of these parameters on required torque and provide insight on how timing may change.

Similar to other undulatory animals moving through water (Gillis, 1998a), the EMG duty cycle decreased in the posterior direction (Fig. 6D) from 0.52±0.07 at 0.3 SVL to 0.40±0.06 at 0.9 SVL. This decrease in duty cycle was small relative to that in some aquatic fish. Studies have shown that for elongate animals with uniform body shape, such as the eel (Gillis, 1998a) and lamprey (Williams et al., 1989; Wardle et al., 1995), EMG duration changes are small along the length of the body compared with those in fish with differing body diameter and shape (Wardle et al., 1995; Altringham and Ellerby, 1999). Also, the duty cycle in sandfish was comparable to that in other animals such as the terrestrial snake [~0.45 in iliocostalis of Nerodia fasciata) (Jayne, 1988)] and aquatic lamprey [~0.5 (Williams et al., 1989)], despite differences in environment. Currently, we do not understand what factors influence this timing and how the environment might play a role. In future studies, we will use models in order to investigate timing and duration as it relates to movement through different environments.

Many studies (Frolich and Biewener, 1992; Biewener and Gillis, 1999; Nishikawa et al., 2007) have found a change in neuromuscular strategy with changing external resistance, suggesting sensory feedback is used to modulate gait during swimming. We argue that the sandfish controls its spatial form to target a template to optimize escape into the granular medium. (1) EMG intensities increased with depth (and thus resistive force), which suggests that the sandfish is sensing changes in external environment. (2) No differences between A/λ at varying depth and ϕ were detected, suggesting that the change in neural strategy allowed the animal to maintain the waveform within a range. (3) RFT and computational models (Maladen et al., 2011b) of an undulatory swimmer in granular media suggest that within this range of A/λ, speed is maximized and mechanical cost of transport is minimized, and these parameters are insensitive to changes in A/λ≈0.2 (Fig. 8). Therefore, targeting this template can yield benefits in speed and energetics as well as sensitivity to changes in kinematics during escape. However, to confirm that this is a general template for undulatory locomotion through sand, kinematics in other sand-swimming animals must be assessed. We predict that sand-adapted subarenaceous animals (defined as animals that primarily live subsurface) with similar body shape should employ a similar waveform to achieve rapid burial. Promising candidates include Sphenops sepsoides and Neoseps reynoldsi.

Unlike the spherical particles used in our experiments, some natural sands have irregular shapes. These irregularities can allow individual particles to interlock, effectively increasing particle–particle friction (Cho et al., 2006). Consequently, volume fractions in natural sands can be lower (Dickinson and Ward, 1994; Cho et al., 2006) than our laboratory LP state; however, the volume fractions we achieve are within the range of volume fractions found in nature. Also, the mechanical properties of spherical glass beads have similar characteristics as natural sands (Goldman and Umbanhowar, 2008): resistance force is insensitive to speed and increases with depth (Maladen et al., 2009). A previous study of sandfish showed that η is independent of particle size over an order of magnitude (Maladen et al., 2011b). Also, simulation of sandfish has shown that grain–grain friction has a small effect on sandfish η while body–grain friction has a larger effect (Maladen et al., 2011a). Increasing particle–particle friction by 50% caused a 14% increase in η while increasing particle–body friction by 50% caused an 84% decrease in η. Thus, we expect that our laboratory granular media elicit similar η and template as those used during swimming in natural desert sand.

The sandfish lizard is a sand-swimmer, displaying both morphological adaptations (Mosauer, 1932; Arnold, 1995; Schleich et al., 1996) and kinematics that allow it to move effectively within granular media (Maladen et al., 2011b). A diversity of animals live on and bury within granular substrates, such as the Namib desert mole (Gasc et al., 1986), sand-swimming shovel-nose snake (Norris and Kavanau, 1966), ocellated skink (Arnold, 1995), caecilian (Summers and O'Reilly, 1997), worm (Dorgan et al., 2005), sand lance (Gidmark et al., 2011) and razor clam (Winter and Hosoi, 2011). The techniques developed for the sandfish (like high-speed subsurface X-ray imaging, subsurface EMG, and substrate preparation using a fluidized bed) could also be used to quantify and compare kinematics for other species and determine whether the neuromechanical principles that govern sandfish movement are similar for other subsurface locomotors.

We found that the muscle activation pattern in the sandfish was correlated with drag forces in granular media, and surprisingly that the general features of muscle activation timing was similar to those of undulatory swimmers in water despite differences in environment. The discovery that muscle activation strategy changes with depth (Fig. 4) while spatial waveform remains similar in sandfish suggests that these kinematics are targeted and maintained using sensory feedback. This finding supports the hypothesis that operating with A/λ≈0.2 is a template for single period sinusoidal undulatory locomotion through sand and essential to attain a high-speed escape with a low cost of transport.

Finally, we note that the sandfish is an attractive animal model on which to test hypotheses and model neuromechanics of locomotion. First, the sandfish swims in granular media over a range of speeds and resistive force conditions while maintaining similar kinematics. Most studies have shown a change in neural activation and kinematics with differing environments (Fang-Yen et al., 2010) and speeds (Gillis, 1998a), which makes their coupled effect harder to interpret. In addition, the sandfish has similar wave amplitude along the length of its body (unlike most undulatory swimmers in fluids) and its shape during swimming can be described well with a sinusoidal wave. The granular environment in which the sandfish moves is also a controllable environment. The forces within the granular media can be tuned by changing surface pressure, compaction of the media, friction between grains and bead density. Weakening the material using airflow (Abate and Durian, 2005; Daniels et al., 2009) can allow the animal to move in a different regime where material inertia becomes important. Such manipulations could provide more insight into activation timing and sensory systems used during swimming.

In future studies, we will use both RFT and numerical simulation (Maladen et al., 2011b) of the sandfish and granular media to determine how internal forces in the animal body relate to the environment. Changes in muscle activation pattern do not always reflect internal force changes and such models can be used to estimate forces during sand-swimming. In addition, our sandfish-inspired robotic model can be used to test neuromechanical control strategies in real-world granular environments (Maladen et al., 2011c). For example, we can change morphology and passive stiffness to test how they would affect neuromechanical phase lags. Also, control systems that modify the robotic model's motor output based on the resistive force can be implemented to optimally change gait with varying conditions. These biologically inspired activation strategies can also give insight into challenges associated with the movement of robots through complex terrain (Choset et al., 2005; Maladen et al., 2011b; Maladen et al., 2011c).

We thank Stephen P. DeWeerth, Paul Umbanhowar, Chen Li and Ryan Maladen for discussions. We thank Mateo Garcia for assistance with side-view sandfish experiments. We thank Thomas J. Burkholder for assistance with sandfish dissection and the Robert E. Guldberg's lab for assistance with a CT scan to identify bone morphology. We would like to thank Elizabeth A. Gozal for providing the EMG analysis code (SpinalMOD). We thank Tonia Hsieh and an anonymous reviewer for advice regarding biological statistics.

 
• A

  wave amplitude

 
• A/λ

  a dimensionless parameter which characterizes the spatial form of sandfish kinematics

 
• CoT

  cost of transport

 
• CP

  closely packed granular medium

 
• EMG

  electromyogram

 
• f

  undulation frequency

 
• I

  EMG intensity – the ratio of the EMG integrated area to the EMG duration

 
• I_rel

  relative EMG intensity – normalized EMG intensity using average intensity during undulation 3 for a single experiment

 
• L_M

  normalized muscle length estimated from the angle between segments

 
• LP

  loosely packed granular medium

 
• RFT

  resistive force theory

 
• STL

  snout-to-tail length

 
• SVL

  snout-to-vent length

 
• ν_c

  speed of the traveling curvature wave

 
• ν_EMG

  average speed of traveling EMG wave

 
• ν_for

  forward speed – speed of a marker in the forward direction of movement of the animal

 
• ν_lat

  lateral speed – speed of a marker in the lateral direction of movement of the animal

 
• ν_seg

  segmental speed – speed of a marker in the direction of movement

 
• ν_w

  speed of lateral displacement

 
• ν_x

  average forward speed

 
• η

  wave efficiency – the ratio of the average forward speed of the animal to the speed of the kinematic traveling wave

 
• θ_d

  angle of descent – angle between horizontal and the best linear fit to the positions of all tracked dorsal markers during the subsurface descent

 
• λ

  wavelength of the sinusoidal wave

 
• ϕ

  volume fraction – the ratio of the solid volume of the grains to the total occupied volume