Locomotor benefits of being a slender and slick sand swimmer

The prediction of locomotor performance in relation to features of body morphology (body shape, skin properties and limb number) is a topic of interest from bio/neuromechanical (Aristotle, 350 bce; Dayton et al., 2005; Gans, 1975; Mosauer, 1932; Russell, 1917; Tytell et al., 2010; Wainwright and Reilly, 1994) and robotic perspectives (Brooks, 1992; Lauder et al., 2007; Maladen et al., 2011a; Maladen et al., 2011c; Pfeifer et al., 2007; Yim et al., 2007). From an evolutionary perspective, it is important to quantitatively address questions associated with putative selective pressures for locomotor apparatus and bauplans in different environments (Lauder, 1981; Lauder and Drucker, 2004; Walker and Westneat, 2000). Such studies are necessarily complicated because body morphology, kinematics and environment interactions are typically strongly coupled.

Much progress in disentangling these factors has been made in the studies of locomotion in fluids. For example, in swimming, using a variety of tools [(e.g. high-speed imaging, particle image velocimetry (Drucker and Lauder, 1999)], many biological studies have been conducted to understand how aspects of morphology [e.g. effects of body shape, limb shape and stiffness (Fish and Battle, 1995; Long and Nipper, 1996; Miklosovic et al., 2004)] and movement mechanics (Blake, 2004; Fish and Lauder, 2006; Liao, 2007; Sfakiotakis et al., 1999) contribute to performance. Such systems have also been the focus of models that allow systematic understanding by linking shape to function, through probing of fluid–structure interactions and how aspects of morphology affect performance. These modeling approaches include mathematical and computational (Alexander, 2003; Anderson and Wendt, 1995; Lighthill, 1971; Liu et al., 1996; McHenry and Jed, 2003; Pedley and Hill, 1999; Tytell and Lauder, 2004; Webb and Weihs, 1986) modeling using the Navier–Stokes equations of fluid dynamics, as well as a more recent modeling approach, using physical robot swimming and flying models (Brooks, 1992; Lauder et al., 2007; Maladen et al., 2011a; McHenry et al., 1995; Yim et al., 2007).

Seemingly complex dry granular environments (which display solid and fluid-like aspects) are also amenable to detailed experiment/theory comparison (Ding et al., 2012; Ding et al., 2013; Maladen et al., 2011b). High-speed X-ray imaging and EMG studies of the biomehanical and neuromechanical aspects of the sand swimming sandfish lizard, Scincus scincus (Linnaeus), a ‘subarenaceous’ animal (Mosauer, 1932) (living within sand and swimming for long distances) native to the Sahara Desert (Fig. 1C) have provided insight into the locomotor strategy and performance of this sand swimmer, and have even revealed principles of swimmers in real fluids (Ding et al., 2013). X-ray imaging of the sandfish lizard (Maladen et al., 2009; Sharpe et al., 2013) revealed that it swims subsurface using single-period anterior-to-posterior traveling body curvature waves while its limbs remain along its body. A multi-particle discrete element method (DEM) simulation of a granular material (composed of 3 mm diameter particles) coupled to a numerical model of the animal captured swimming kinematics and predicted that sandfish swimming takes place within a self-generated localized ‘frictional fluid’ (Ding et al., 2012). However, because of the high computational time, DEM could not be used to study locomotion in smaller particle natural substrates.

A granular resistive force theory (RFT) was developed (Ding et al., 2012; Maladen et al., 2009; Maladen et al., 2011b) to gain insight into the principles of sand swimming in more natural environments. In RFT, a body is partitioned into infinitesimal segments, each experiencing thrust and drag related to the orientation and environmental interaction. RFT assumes forces on a particular segment are independent of forces from neighboring

segments. Net forward thrust is determined by integrating forces over all body segments. With empirically measured force relationships (Maladen et al., 2009) (the granular equivalent of Stokes' laws), the granular RFT allowed us to make accurate statements about sandfish swimming kinematics and performance. The model revealed that its body waveform (i.e. curvature) during swimming enables the sandfish to operate at maximal undulation efficiency (Maladen et al., 2011b) while expending minimal work (Sharpe et al., 2013). This indicated that the waveform selected was a target of control, or a ‘template’ (Full and Koditschek, 1999), which muscle activation experimentation confirmed (Sharpe et al., 2013). Furthermore, RFT was successfully used to make predictions of neuromechanical activation patterns during swimming (Ding et al., 2013). This was a consequence of the simplicity of the sandfish shape and body wave, the lack of inertial and non-local effects, and the dominance of environmental forces relative to internal forces.

Compared with other subarenaceous animals that move primarily in dry sand, the sandfish has a relatively small ‘elongation ratio’ of L/w=7, which we define as the ratio of L, the effective body length (effective length accounts for tail tapering, see Materials and methods), to w, the average body width; see examples in Table 1. Other dry-sand subarenaceous squamates exhibit larger elongation, with Neelaps bimaculatus having L/w≈50. Elongation of the body among sand and soil burrowers has been hypothesized to enhance subsurface locomotion performance and possibly energy efficiency (Gans, 1975; Lee, 1998).

In this study, we used the tools described above to test this hypothesis, in particular to study the effects of elongation ratio on sand swimming performance. We did so by performing the first quantitative exploration of the kinematics of locomotion of the shovel-nosed snake, Chionactis occipitalis (Hallowell) (Fig. 1A), a desert-dwelling snake whose L/w≈30 is typical of subarenaceous swimmers. We then used RFT to demonstrate the advantages gained by being long and slender particularly in terms of the efficiency of movement and mechanical cost of transport. RFT allows us to decouple and independently examine the effects of kinematics, body shape and skin–grain friction on sand–swimming performance, all of which could impact locomotion. Therefore, RFT enables us to investigate the extent to which body elongation plays a role. While the sandfish and shovel-nosed snake are evolutionarily constrained by different bauplans (i.e. have divergent body forms and skin properties), our results demonstrate that they have achieved optimized (but not convergent) movement patterns for subarenaceous locomotion.

Morphological measurements were made on four anesthetized shovel-nosed snakes and 11 anesthetized sandfish. Sandfish have 26 presacral (trunk) vertebrae and well-developed limbs (Sharpe et al., 2013), whereas shovel-nosed snakes are limbless and have a variable number of precaudal vertebrae (range 165–175), with our animals representing a population near Yuma, AZ, USA (Cross, 1979). Importantly, with respect to the independent convergence on an elongate and limbless morphology and burrowing ecology highlighted by Wiens et al. (Wiens et al., 2006), our comparison here involves two species that swim through sand, but do so with decidedly different morphologies.

An effective body length was defined for both animals to account for tail tapering (see Materials and methods). L/w for the snake was ~5 times larger than that for the sandfish (33.3±2.5 and 7.0±0.4, respectively, F_1,13=1348, P<0.01, Table 2). The winding number (W) is defined as the number of coils an animal can form along its body when maximally bent and is equivalent to W=(κ_max/L)/2π, where κ_max is the maximum bending curvature. We measured W on anesthesized animals (three snakes and three sandfish) to bound the maximum number of waves attainable for a given curvature (Fig. 1B,D). The average winding number was larger for the snakes than for the sandfish, W=6.2±0.8 and W=1.4±0.06, respectively (F_1,4=114, P<0.01). Snake body–particle friction (see Materials and methods for experimental details) was ~50% lower than sandfish body–particle friction on the ventral surface (forward direction: snake μ_s=0.11±0.02, sandfish μ_s=0.19±0.02, F_1,4=29.6, P<0.01; backward direction: snake μ_s=0.14±0.02, sandfish μ_s=0.28±0.03. F_1,4.1=61.9, P<0.01; Tables 2 and 3). However, body friction was similar for the dorsal surface of the snake and the sandfish (forward direction: snake μ_s=0.16±0.03, sandfish μ_s=0.18±0.03, F_1,4=0.6, P=0.48); backward direction: snake μ_s=0.24±0.03, sandfish μ_s=0.23±0.02, F_1,4.1=0.04, P=0.85; Table 3).

We used a fluidized bed (see Materials and methods and Fig. 7A) to prepare the granular medium into a loosely or closely packed state (LP or CP). [Note that as kinematics were statistically similar (see supplementary material Table S1) during swimming in these different states, we will report the combined LP and CP averages.] After a gate separating the snake from the granular media was lifted, snakes moved from the rigid surface of the holding container and buried into the prepared sand; during the initial stages of burial, part of the body was in contact with the holding container (Fig. 2A). Snakes buried themselves using an approximately serpenoid waveform with lateral body undulations that propagated from head to tail (Fig. 2B, Fig. 3A; supplementary material Movie 1). During subsurface sand swimming, X-ray imaging revealed that a wave of curvature propagated down the body at 3.7±1.3 cm s⁻¹, 10 times slower than sandfish swimming in which waves of curvature propagated at 39.5±15.9 cm s⁻¹. We estimated the angle of entry of the snakes (the average downward angle of the body relative to the horizontal) using the known arc length between markers and the measured projected distance between markers acquired in the top-view X-ray. We found that snakes moved into media at an entry angle of 23±10 deg. This value was similar to the angle of entry found for the sandfish (Sharpe et al., 2013).

On initial observation, the shovel-nosed snake appeared to employ a different locomotion strategy from the sandfish – that of tube following, similar to the nematode worm Caenorhabditis elegans on an agar surface in which a groove is formed and the material does not yield around the body (Berri et al., 2009). However, when comparing the animal's midline position over a few undulation cycles (Fig. 3A), it was evident the snake did not move in a perfect tube. Previously, we quantified such ‘slip’ using a measure called the wave efficiency, the ratio of the forward speed to the wave speed (Maladen et al., 2009). This measure assumes that wave speed is constant, the animal's body curvature does not change, and the animal moves forward, on average, along a linear trajectory. Because the snakes often turned during sand swimming such that the body shape was more complicated than a constant amplitude (and direction) wave, we defined a new parameter called the local slip angle, β_s, which is a measurement of the angle between the direction of movement of a body segment and the body segment orientation. We measured this parameter at every marker along the body during subsurface movement using the relationship , where is the unit velocity vector and is the tangent vector at a body marker location (m) and time (t). Physically, this quantitatively measures the local deviation from ‘tube’ movement.

The slip measurements depend on the accuracy of the vector calculations. To determine body midline positions and tangent vector directions in the presence of relatively large detector noise (occlusions from above-surface chamber walls, and movement of the granular media), we developed a marker detection and tracking strategy that provided the trajectories of multiple markers to sub-pixel accuracy. We used a recursive Bayesian approach to marker trajectory estimation with linear models and Gaussian noise, leading to a combined Kuhn–Munkres algorithm and Kalman filtering formulation. The inclusion of a B-spline model for predicting the subject's movement provided additional robustness to image artifacts and unobservable markers. For details, see Materials and methods and supplementary material Movie 1 (M.S., S.S.S., D.I.G. and P.A.V., in preparation).

We define the spatial and temporal average of the slip angle as ; the average is computed over all markers visible during a trial. The markers, m, over which β_s was averaged varied depending on marker visibility, but we always excluded the first and last marker, which were more susceptible to error due to the B-spline fitting error at the end points. is inversely related to the wave efficiency for a sinusoidal straight swimmer, but is more general because it can be calculated for animals that move along any arbitrary path, change undulation frequency, or change shape during movement.

was measured for snake swimming (three snakes, 18 trials) and sandfish swimming (six sandfish, 56 trials). Angular differences between velocity and tangent vectors at the tracked marker positions (Fig. 3C,D) varied along the body and during time for both animals but were larger during sandfish subsurface locomotion compared with snake locomotion, resulting in high local slip (β_s) measurements (see Fig. 3E,F, for β_s during a representative snake and sandfish trial, respectively). Thus, the average slip angle for snake sand swimming was small (=6.4±1.4) compared with sandfish sand swimming (=21.3±4.1, ANOVA, F_2,73=280, P<0.01).

Previous work suggested that the sandfish wave efficiency (and therefore slip) depended on kinematic variables (Maladen et al., 2011b) such as maximum curvature and undulation number. Therefore, we next examined differences in kinematics between the snake and sandfish by measuring the relative curvature and number of waves along the body during sand swimming for each. We measured the maximum local curvature, κ_max, using a best fit circle to each bend and measured the half-arc length, λ_s/2, determined as the length between locations of zero curvature, for each bend (see Fig. 4A); here s denotes arc length. The relative curvature, κλ_s, determined as 2κλ_s/2 of each bend along the body, was calculated for each snake and sandfish trial at an instant in time when most (if not all) of the animal was within the field of view and for submerged body portions only. κλ_s did not vary with media compaction for either animal. Snakes had a broader distribution of κλ_s during swimming compared with sandfish, where κλ_s of each bend along the body was considered an independent measurement (Fig. 4C). The mean κλ_s () for a single trial was calculated by averaging κλ_s over all bends visible along the body.

On average, snake trials had a slightly but significantly lower than sandfish trials (=6.1±1.3 for snakes, and =7.4±1.0 for sandfish; ANOVA, F_2,73=20.3, P<0.01; Table 2). Because the snake was more inclined to turn during movement, we compared the distribution of κλ_s for straight swimming trials, in which the snake was moving on average along a straight trajectory, and for all swimming trials, in which the snake may turn or move straight (see supplementary material Fig. S2). Both sets of trials showed a broad distribution of κλ_s, indicating that snakes employ small and large κλ_s even during straight undulation cycles. Furthermore, we found the curvature along each bend decreased as the snake moved farther (and therefore deeper) into the media (ANOVA, F_7,105=3.3, P<0.01 for all trials considered), but this trend was more pronounced when considering straight swimming alone. For straight swimming, the maximum κλ_s decreased as the snake moved deeper into the medium (see Materials and methods and supplementary material Fig. S2).

The snake also swam with a higher number of waves along its body; we refer to this quantity as the undulation number, ξ, determined from counting the waves on the animal's body at an instant in time when all or most of the body was within the field of view (see Materials and methods for more details). The mean undulation number across all trials, , was larger for the snake (=3.5±0.7) than for the sandfish (=1±0.1, ANOVA, F_2,73=695, P<0.01) and did not statistically vary with granular media compaction for either animal (Table 2). As both wave curvature and wavelength changes along the body, we explored local characteristics by calculating a related measure called local undulation number, defined as the ratio of effective body length to twice the arc length of each body bend (or half wave), L/λ_s; if λ_s is constant along the body, then L/λ_s=ξ. Similar to ξ, the distribution of L/λ_s for the snake was wider and the average was larger than that of the sandfish (Fig. 4D).

The low slip locomotion of the snake appears ‘efficient’ compared with that of the sandfish. However, it is unclear whether these kinematic differences are due to lower skin friction, the larger number of waves propagated along the body during swimming, the snake's more slender body, or some combination of these factors. Inspired by our observation that the slip is finite and thus material yields in a small region around the body, we next applied the granular RFT (discussed in Materials and methods) to the snake to tease apart the effect of these factors on subarenaceous movement. In the model, we assumed that the animal targeted a particular kinematic traveling wave pattern, thereby applying a neural control strategy to bend its body into a desired shape; we will return to this point later in the paper. To calculate slip on such a deforming shape, we generated waves of constant curvature, wavelength and undulation number, and simulated horizontal plane swimming at a depth of 7.6 cm (see Materials and methods). We modeled swimming for two different body morphologies: one with elongation ratio L/w=7.1 (similar to the sandfish) and another with L/w=33.5 (similar to the snake), but with the same coefficient of body–particle friction, μ_s=0.17. Note that both turning and straight (i.e. average midline follows a straight path) snake experimental swimming trials were considered in our analyses. Analysis of all trials was supported by theory, which revealed that simulation of a swimmer moving in a circle with a given average κλ_s yielded similar β_s values to those for straight swimming.

The theory predicted that increasing curvature resulted in decreased slip, (Fig. 5). Using a larger undulation number () also decreased but yielded diminishing returns for >2. Despite the reduction in slip due to the long, slender body, the RFT calculations did not match the observed slip values for the snake. However, because the snake's ventral skin–grain friction was approximately half that of the sandfish, we next decreased friction in the RFT simulations by reducing tangential drag forces by 50% (see supplementary material Fig. S1) (Fig. 5, dashed curves). Theory predicted a greater than 50% decrease in , and thus fell within the experimentally observed slip range for the snake when considering its kinematics.

Our previous work demonstrated that the sandfish targeted optimal kinematics to maximize undulatory efficiency (η_u, forward distance moved per cycle normalized by λ_s) and minimize the mechanical cost of transport (CoT, sum of the power to external surroundings, normalized by body weight and speed, neglecting internal forces) (Sharpe et al., 2013). These results suggested that intermediate curvature waves can be a control target for optimal sand swimming and could therefore be a ‘template’ for subarenaeous locomotion. In accord with the optimal swimming hypothesis, we found that the snakes also operated within a and range that is optimal (Fig. 6A,B). As with previous work exploring sandfish swimming efficiency (Maladen et al., 2009), there is a slight difference between undulation efficiency in LP and CP media. There is a larger difference in CoT; swimming in LP requires approximately half the power. However, the minimum CoT occurs when using similar for both CP and LP simulations. Accordingly, we found statistically similar CP and LP experimental kinematics for both sandfish and snakes (see supplementary material Table S1 for separated CP and LP values). We note that simplifications necessary in modeling and/or physiological pressures could account for the slight difference between experimentally observed kinematics and the model-predicted optimal kinematics.

Experiments revealed that the snake swam in a frictional fluid with a lower average slip angle than the sandfish. Granular RFT captures this performance differential and revealed at least two locomotor benefits of having a long, slender body. (1) For a given curvature and undulation number, a larger elongation ratio yields a slightly smaller slip. This suggests that when using the same kinematics, a more slender and/or longer body will yield lower slip locomotion. However, the change in slip is small and cannot account for the large decrease observed between the sandfish and snake. (2) A larger winding number (a consequence of possessing a larger L/w and a greater number of vertebral segments) allows an undulatory swimmer to increase its undulation number while maintaining a high average curvature; therefore, animals with a larger elongation ratio (and an increased number of vertebrae) can utilize a greater range of curvatures and undulation numbers. In contrast, an animal with a low winding number may have to decrease curvature to achieve a higher undulation number and vice versa. For example, for the sandfish's body shape (L/w=7.1) and winding number (W=1.4), if =2 the swimmer is constrained to operate with <4.3 (Fig. 5, solid cyan curve). Therefore, it is more advantageous for a swimmer with a low elongation ratio to utilize a smaller number of waves along its body, which enables high curvature swimming kinematics. However, for the slender snake, this compromise is not necessary (Fig. 5, see dotted-dashed cyan curve).

Our results also imply that the longer body and lower friction allows the snake to achieve a higher undulation efficiency and lower mechanical CoT than possible for a sandfish given its bauplan (Fig. 6). We postulate that the bauplan exemplified by the shovel-nosed snake may also approach that of an optimal sand swimmer. While a longer body with a greater number of vertebrae would allow the animal to utilize a higher undulation number or relative curvature, RFT calculations show diminishing returns in undulation efficiency (Fig. 5 and Fig. 6A). In addition, a longer body would also incur greater frictional drag, and thereby increase CoT. Further systematic exploration of the RFT model is needed to test this idea.

Our results suggest how low skin friction offers a biomechanical advantage for sand swimming and support prior statements that smooth skin is an adaptation that is characteristic of subterranean animals (Mosauer, 1932). However, more quantitative studies that measure skin friction using the same technique and on a multitude of species, both fossorial and non-fossorial, are needed to support this claim. And while RFT predicts that lower friction would improve performance, the shovel-nosed snake skin friction (μ_s≈0.11 on the ventral surface) may be one of the lowest found for squamates in nature (Berthé et al., 2009; Gray and Lissmann, 1950; Hu et al., 2009) and additional reductions could be morphologically constrained.

Somatic elongation has evolved independently many times among squamate reptiles, involving all snakes and multiple examples among the ‘lizards’ (i.e. non-snake squamates) (reviewed by Wiens and Slingluff, 2001). Elongation has been achieved differentially by elongation of the trunk and/or the tail, and typically includes reduction or loss of the limbs. Burrowing ecology has been associated with elongation in squamates, but many elongate taxa do not burrow. The burrowing bauplan typically includes elongation of the trunk with a relatively short tail (Wiens et al., 2006), and this association holds for the shape of the shovel-nosed snake. Thus, the convergence of a long trunk (via an increased number of precaudal vertebrae), relatively short tail and reduction/loss of limbs in association with burrowing is a well-established repeated theme in the evolution of squamates (Brandley et al., 2008; Sites et al., 2011). Our comparison here focused on two species that have converged ecologically on a subarenaceous burrowing lifestyle, yet unlike typical evolutionary comparisons (e.g. Wiens et al., 2006) the sandfish and shovel-nosed snake are not convergent in morphology. The shovel-nosed snake is a typical serpent, albeit with a relatively short tail, while the sandfish is a typical limbed lizard that is phylogenetically nested in a clade that does not include elongation or limblessness (Wiens et al., 2006). Furthermore, the sandfish lizard uses its limbs during initial burial to achieve rapid submergence and during foraging activities. Our comparisons reveal evident patterns of performance optimization and constraint based on two remarkably different squamate body forms. Furthermore, this study reveals the functional benefits of elongation on burrowing locomotion, providing the first quantitative support for the hypotheses made by Gans (Gans, 1975) that elongation leads to enhanced subsurface locomotor performance and decreases energy expenditure, and provides a functional link between correlations noted between burrowers and adaptive body features (Lee, 1998).

We used X-ray imaging to reveal the kinematics of sand swimming of the shovel-nosed snake, and then used a granular RFT to quantitatively determine how the snake's body morphology and skin friction conspire to achieve high performance. Our model reveals that a long, slender, slick (and limbless) (Lee, 1998; Wilson, 2012) body leads to high swimming performance in subarenaceous animals – the snake's long body and high winding number allow it to operate over a wide range of curvatures and wave numbers. Both subarenaceous animals in this study swim with a body wave curvature that maximizes undulatory efficiency and minimizes the mechanical CoT for their particular bauplans. This provides evidence that both animals target these kinematics, supporting the idea that control for such waveforms is a ‘template’ for sand swimming.

More broadly, sand swimming continues to be an excellent system for making detailed comparison between experiment and theory in a locomoting biological system. In this non-inertial, friction-dominated environment, internal body forces are small compared with environmental forces (Ding et al., 2013), and shape changes determine locomotor efficacy. Combined with the relative simplicity of the granular frictional fluid environment, the system allows the teasing apart of relative contributions of shape, friction, wave properties and even aspects of neuromechanical control (Ding et al., 2013). Such tools might also prove useful to address questions in functional morphology and ecologically driven adaptation (Garland and Losos, 1994) in dry granular media (e.g. cranial modifications, smooth skin and miniaturization) that occur in a diversity of burrowing reptiles (Gans, 1975; Lee, 1998; Mosauer, 1932; Rieppel, 1988; Wiens and Slingluff, 2001; Wilson, 2012).

Shovel-nosed snakes are found in the Mojave and Sonoran deserts, typically near dry sand dunes (Norris and Kavanau, 1966). The shovel-nosed snake moves only short distances on the surface of sand, in comparison with other snakes, and readily buries into granular media (Mosauer, 1932; Mosauer, 1933; Mosauer, 1936; Norris and Kavanau, 1966) to depths of 60 cm beneath the ground (Cowles, 1941). Chionactis occipitalis were collected in accordance with a scientific collection permit (no. SP591773) approved by the Arizona Game and Fish Department. All snake experiments were conducted under the Georgia Institute of Technology IACUC protocol A11066.

Sandfish lizards, Scincus scincus, were purchased from commercial vendors (LLL Reptile, Vista, CA, USA and Ocean Pro Aquatics, Chino Hills, CA, USA). Data sets were analyzed from previous work (see Sharpe et al., 2013). All experimental procedures were conducted in accordance with the Georgia Institute of Technology IACUC protocol number (A08012) and Radiation Safety protocol (X-272).

Morphological measurements were made on four anesthetized snakes and 11 anesthetized sandfish (Table 4). Both animals were sedated by 5% isoflurane gas induction and then anesthetized with an intramuscular injection of ketamine hydrochloride (200 mg/kg body mass).We report average masses and lengths for animals that were measured more than once over the course of the study.

SVL, STL, w and mass were recorded for each animal. Snakes had a consistent body width, w≈1 cm, and long bodies with an average STL of 36.6±2.8 cm (see Table 4). A noticeable reduction in body diameter occurred between 0.85 and 0.9 STL. Therefore, we defined an effective snake length (L) to account for this tail tapering as L=0.87STL. On average L=31.8±2.4 cm for snakes. The sandfish width was w≈1.5 cm, the STL was on average 14.1±0.9 cm, and the SVL was 8.9±0.5. We accounted for tail tapering in sandfish by defining L as 1.2SVL, determined from morphological observations, which was on average 10.6±0.6 cm. Snake and sandfish masses were similar (17.5±3.4 and 17.6±2.8 g, respectively). The elongation ratio, L/w, was calculated using the effective length and body width measurements. L/w for the snake was ~5 times larger than that for the sandfish (33.3±2.5 and 7.0±0.4, respectively).

Morphological measurements were also made on awake animals used in kinematic experiments on the day of data collection. Because awake animals tended to bend their bodies and move during measurements, the measurements were subject to more error in accuracy. However, average L/w values were similar for these animals to those made on the anesthetized animals (35.6±3.4 for the snakes, N=3 animals, and L/w=7.3±0.2 for sandfish, N=6 animals).

Winding number (W) was measured by bending the bodies of three anesthetized snakes and three anesthetized sandfish lizards. Bodies were maximally bent to determine the total number of coils that could be formed along the body from the snout tip to the effective body length. W was used to determine the maximum relative curvature that the animal could attain given a serpentine waveform and fixed number of waves along the body.

Body–particle friction measurements were made using 0.3 mm glass particles (the experimental media) following techniques described elsewhere (Hu et al., 2009) (Fig. 7B). Animals were placed on an inclined plane covered with a monolayer of particles such that either their ventral or their dorsal surface was in contact with the particles. Their head was either oriented downward (causing ‘forward’ sliding) or upward (causing ‘backward’ sliding) relative to the sloping angle. In total, four orientations were tested. The edge of the plane was lifted until an inclination angle (θ_slide) was reached that caused animals to slide (Fig. 7B). The static coefficient of friction, μ_s, was calculated from the relationship μ_s=tan(θ_slide). An anesthetized animal was placed on the monolayer of particles adhered to a board. One edge of the board was lifted and θ_slide at which animal slid was recorded. Four friction measurements were made in each orientation (unless otherwise noted in Table 3) on three snakes and three sandfish lizards.

LP and CP compactions were achieved using a fluidized bed and shaking apparatus similar to prior experimentation (Fig. 7A) (Gravish et al., 2010; Li et al., 2009; Sharpe et al., 2013). The bed (22.9×40.6 cm²) was filled with glass particles (diameter=0.27±0.04 mm and density=2.5 g cm⁻³) to a height of 10 cm; such particles approximate fine desert sand. The same bed and granular particles were used for both snakes and sandfish kinematic trials, but slightly different shaking apparatuses were used to generate CP states. Consequently, a 2% difference in CP volume fraction (ϕ, measured as the ratio of the volume of the granular medium to the volume of its occupied space) occurred between snake and sandfish experiments, where snake CP trials were at ϕ=0.617±0.012 and sandfish CP trials were at ϕ=0.635±0.013. A LP compaction of ϕ=0.58 was achieved for both sandfish and snake experiments. Because compaction did not influence kinematics for either animal (see supplementary material Table S1), the difference in CP states should not affect the results.

Snake experiments took place at an average ambient temperature of 25°C and average humidity level of 41%. Lights heated the granular substrate for ~1 h prior to experimentation and the temperature of the sand was measured prior to each trial. The average sand temperature for all trials was 32.5±2.5°C. The bed was periodically fluidized, which helped homogenize the distribution of temperatures within the bed, and lights were turned off when necessary to keep a constant sand temperature. A similar setup was used for the sandfish lizard experiments (see Sharpe et al., 2013).

Previous sandfish locomotion studies showed that sandfish kinematics were similar in the two preparations (Maladen et al., 2009; Sharpe et al., 2013). Similar snake kinematics were obtained in LP and CP preparations; therefore, the results reported represent the combined data for the two preparations unless otherwise specified (see supplementary material Table S1 for separated kinematic results for LP and CP trials).

Subsurface kinematics were recorded for three of the snakes (snake nos 1, 3 and 4; see Table 4) and each were tested in three CP and three LP granular media preparations, with a total of 18 snake trials considered in the analysis.

Kinematic sandfish data were utilized from previous studies (Maladen et al., 2009; Sharpe et al., 2013). Fifty-six undulations were considered from different trials taken across six animals for comparison with snake data. All experimental procedures were conducted in accordance with the Georgia Institute of Technology IACUC protocol nos A11066 and A08012 and radiation safety protocol no. X-272.

Snakes were placed in a container resting above the surface of the prepared sand (Fig. 7A). A gate, separating the snake from the granular media, was lifted and the bottom-rear section of the container was also lifted via a string connection, causing the snake's orientation to angle downward toward the media. Transparent walls resting on the surface of the sand outside the holding container ensured burial at a desired location away from container walls and within the X-ray field of view. Above-surface visible light video was recorded at 63 frames s⁻¹ using a high-speed camera (AOS Technologies AG X-PRI). Subsurface snake kinematics were recorded using an X-ray system (OEC 9000) coupled to a high-speed camera (Fastcam 1024 PCI, Photron). Top-view (i.e. dorsal view) subsurface snake kinematics were recorded at 60 frames s⁻¹. The X-ray system was set to an energy of 85 kV and 20 mA. Because of the apparatus position and X-ray image intensifier size, the field of view within the container was limited to a 20 cm diameter circle; all portions of the snake were not always visible in X-ray images. To enhance contrast, small markers (~1 mm² each, mass ~0.01 g) made from lead were bonded to the midline of the snake using cyanoacrylate glue. These markers were placed at 1 or 2 cm increments down the length of the body. Markers were tracked using custom-written MATLAB code (M.S., S.S.S., D.I.G. and P.A.V., in preparation). If a marker did not exist at every 1 cm location, that location was manually tracked over the course of the trial. A third-order Basis spline (B-spline) was fitted to the snake's body at every frame for all markers visible (see below for more details).

Kinematics were characterized during subsurface movement. The number of undulations along the body, ξ, was determined by counting the waves along the body (to the nearest 0.1 wave) for all animals during a time when most, if not all, of the body was within the field of view. When the full length of the snake was not within the field of view or subsurface, ξ was calculated as ξ_inside+ξ_outside. ξ_inside was directly counted, while ξ_outside was estimated by using the known length of the snake outside of the field of view and the average wave arc length (λ_s). The local undulation number was also characterized and defined as the ratio of the animal's effective length to the local wavelength, L/λ_s. λ_s was found for each bend using X-ray images during sand swimming and measuring the half-wave arc length, λ_s/2, between two adjacent locations of zero curvature and multiplying by 2 (Fig. 4). λ_s was then scaled by the estimated snake entry angle to compare with the effective body length, L. Relative local curvature (κλ_s) was used for a non-dimensional comparison with the sandfish data and calculated as 2κλ_s/2, where κ was determined for each bend using a best fit circle near the maximum. To calculate slip, tangent vectors, , at varying body locations were found by using the tangent to a B-spline that was fitted to the midline position of the snake at every frame (see below). Unit velocity vectors, , were spatially filtered to remove tracking noise. B-splines were also fitted to the sandfish midline during previously recorded trials. Only subsurface kinematics were considered for analyses, where ‘subsurface’ was defined as when, at minimum, the cloacal opening was covered by granular media. Identical methods to those described for the snake were used to estimate ξ, λ_s, κ and for the sandfish.

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.

The granular RFT [based on the model introduced in Maladen et al. (Maladen et al., 2009) and illustrated in Fig. 8] used a simplified model of the snake and sandfish kinematics such that κλ_s values for all bends along the modeled undulatory swimmer's body were the same and no turning occurred. The modeled animals had a uniform body with a circular cross-section and blunt head, and propagated a serpentine wave from head to tail. Bending stiffness measurements in Ding et al. (Ding et al., 2013) demonstrated that for the sandfish, internal forces were small compared with forces required to move through the medium. Therefore, we also used this assumption for the snake.

To numerically solve the force and torque balance on the swimmer, an iterative approach was employed to determine the infinitesimal translation and rotation of the swimmer from an infinitesimal undulation.The force and torque balances provide three equations, which are used to solve for the two in-plane velocities and the rotational velocity. To facilitate the numerics, we use the approximation , and iteratively calculate the motion from the force and torque balances with decreasing w.

To account for the lower friction coefficient of the sand snake skin, the tangential forces of the model were reduced accordingly (see supplementary material Fig. S1). The simulation tracked the instantaneous velocity of each segment, which is used in the slip angle calculations, along with the net displacement and expended work to determine the undulatory efficiency, η_u and CoT. We estimate the resulting numerical error in displacement and work for our results to be no more than 2%.

The number of basis functions defined for each snake was equal to the number of markers along the body (between 26 and 35 for the snake and between 9 and 10 for the sandfish); hence, we were able to fully capture the subjects wide range of motion. These markers were equally spaced along the body at ≈1 cm increments for the snake and 0.1 SVL increments for the sandfish. We used 3rd order B-spline basis functions to have smooth and continuous 1st and 2nd order derivatives along the body. A regularization term λ_b of 0.001 was used for all snake and sandfish data.

Numerical discrete element simulation was used to determine the effect of body-particle friction on parallel and perpendicular forces acting on small intruder moving through an experimentally validated model of 3 mm diameter closely packed glass particles (see Maladen et al., 2011b). A 4 cm-long, 1.6 cm-square rod was dragged such that its long axis was oriented at varying angles (ϕ) with respect to the velocity direction. Only the forces along the length of the rod were considered in analysis; forces on the end-caps were ignored. Body-particle frictions of 0.27 and 0.13 were simulated. As expected, halving the body-particle friction resulted in similar perpendicular forces while parallel forces decreased by 50% (supplementary material Fig. S1).

We thank Yang Ding for performing the drag simulations. Thanks to Perrin Schiebel for providing photos of the Chionactis occipitalis. Thanks to Kevin and April Young for assistance with collection and care of the Chionactis occipitalis.