The non-desert specialist snakes from the collection at Zoo Atlanta were tested in a 2 × 1 m² air-fluidized bed filled with sand collected in Yuma County, Arizona, USA. The snakes’ movement was captured at 30 frames per second (fps) using an overhead camera (data collected for Marvi et al., 2014). The desert-specialists were tested in an air-fluidized trackway of area 152 × 53 cm² filled with 297 ± 40 μm particles (Potters Industries spherical glass Ballotini beads). C. occipitalis moved quickly relative to the non-sand-specialists so its kinematics were captured at 250 fps. In all experiments, a blower attached to the bottom cavity of the fluidized bed was turned on and air flow increased until the GM was in a fluid-like state. The air flow was then decreased until the GM settled in a loose packed state (as described in Maladen et al., 2009). Air flow was always off during experiments and the penetration depth of the snakes (O(mm)) was much less than the total depth of the GM (O(cm)).

We used custom MATLAB code to digitize the midlines of the snakes; image processing functions identified the non-sand-specialist and the method described in Sharpe et al., 2015 tracked the sand-specialist’s banded black markings (Schiebel et al., 2020). A cubic spline fit to the tracked data evenly divided the snake midlines into 100 measurements over the entire body (non-sand-specialists) or from neck to vent (C. occipitalis).

The use of alternating left and right bends, typical of lateral undulation in snakes, on the body was ubiquitous, although specifics of the waveform such as amplitude and number of bends on the body varied (Figure 1). Among the non-specialists, some snakes used uniform, periodically repeating bends (e.g. Figure 1 (b,e,f)) like C. occipitalis (Figure 1 (c,h,i)). Others used bends of varying amplitude and wavelength along the body (Figure 1(g)).

We used the tangent angle, $\theta$, the angle between the local body tangent vector and the average direction of motion of the animal (Figure 2 (a,left)), to represent body posture at each instant in time. Consistent with previous work (Schiebel et al., 2019), we found C. occipitalis kinematics were well-described using a serpenoid curve (Figure 2—figure supplement 1; Hirose, 1993),

Figure 2 with 2 supplements see all

Download asset Open asset

The speed that the wave propagates down the body, $v_{seg}$, determines how quickly the shape changes in time, $t$. At any moment, the shape of the snake as a function of arclength along the body, $s$, is determined by the amplitude, ${\theta }_m$, also referred to as the attack angle, and spatial frequency, or number of complete waves on the body, $\xi$ (Figure 2(a), right). For an infinitesimally thin curve with infinite resolution these parameters are independent. In practice, the musculoskeletal system and width of the body limits both how quickly the midline can change curvature and how sharply bent the body can be, although sufficiently far away from these limits ${\theta }_m$ and $\xi$ remain independent. While there may be other, not purely anatomical, factors which connect these variables, in this work we chose to measure their values and study the connection to performance as mediated by the GM. Future study could explore the neuromechanical origin of the shapes and its connection to the differences observed between different species.

Inspired by this two-parameter description, we characterized snake waveforms using ${\theta }_m$ and $\xi$. We calculated attack angle and wavenumber at each frame in a trial (see Materials and methods for details) and took the mean over all frames to obtain ${\theta }_m$ and $\xi$ representative of the average snake posture. Plotting the results for each of the nine C. occipitalis in the $\xi$ versus ${\theta }_m$ parameter space revealed that the animals used a limited subset of the shapes they were anatomically capable of adopting (Figure 2(b)). There was a slight downward trend in the relationship between $\xi$ and ${\theta }_m$ when considering both each individual and each trial (Figure 2(b), inset). Given the weak relationship between $\xi$ and ${\theta }_m$, we combined all individual measurements and examined the impact of the average C. occipitalis wave, taken as the mean and range of the averages shown (Figure 2(c), black cross).

Compared to C. occipitalis, and consistent with observation (Figure 1), there was greater variety among the wave parameters measured on the non-sand-specialist species (Figure 2(c)). However, as in C. occipitalis, the snakes did not fill the space of anatomically possible shapes. The various combinations of waveforms and performance suggested the forces acting on the snakes depended, at least in part, on the shape used, while the limited space of waveforms alluded to the presence of internal and/or external constraints.

We observed deformations in the granular substrate caused by the motion of the animals (Figure 1(a–c)). This indicated that the shape changes carried out by the animals yielded the media, a situation which necessarily resulted in reaction forces acting on the animal from the GM. We thus explored the connection between body segment motion and the consequent forces by carrying out granular drag experiments using a partially submerged plate.

The snake speeds were variable ($v_{seg}$ from ≈35 to 95 cms^-1, Figure 4—figure supplement 1) and the intrusion depth of the snakes’ trunk into the GM ranged from 0 (no intrusion, occurring at the apexes of the wave which the snake lifted off of the surface) to ≈ 5 mm (Figure 4—figure supplement 2). Previous studies indicated that granular drag stress depends on both intruder speed (Percier et al., 2011) and depth (Albert et al., 1999; Marvi et al., 2014).

Commensurate with those studies, we fixed ${\beta }_d=90$ deg and $z=8$ mm and found normal stress quadratically increased as $v_d$ increased from 1 mm s^-1 to 750 mm s^-1 (the limit of the robot arm capability, Figure 4(a)). Similarly, for ${\beta }_d=90$ deg and $v_d=10$ mm s^-1, normal stress linearly increased with $z$ from 4 mm, the shallowest depth where force could be resolved, to 40 mm, where the plate was fully submerged with the top edge 10 mm below the surface (Figure 4(b)).

Figure 4 with 2 supplements see all

Download asset Open asset

For snakes moving slowly on non-deformable surfaces, body inertia was small compared to the frictional forces between the belly and the surface (Hu et al., 2009). While C. occipitalis was moving quickly, we observed, in line with earlier reports (Gray, 1946), that forward motion would immediately cease when the animals stopped propagating the wave. This phenomenon is observed in swimmers at low-Re where the dominant resistive forces of the surroundings inhibit gliding.

The animal behavior indicated that friction between the body and the GM and local, dissipative interactions within the GM were the dominant forces in the system. Additionally, once the drag experiment was in the steady state, we did not observe any time-dependent structures in the GM. Thus, we used the ansatz that stress on the plate could be modelled as independent, time-invariant static stress and material inertia terms,

The $\rho v_d^2$ term in the formulation is the momentum transfer requirement for pushing the granular material ahead of the plate element, and the static ($v_d\to 0$) stress term, ${\sigma }_o$, is the stress required to initiate yielding of the media. $\rho$ is the density of the material, which we estimated using the density of glass, ${\rho }_{glass}$ and the packing fraction $\varphi =0.64$ induced by motion of the intruder, ${\displaystyle \rho =\varphi {\rho }_{glass}=0.64\times 2500{\mathrm{k}\mathrm{g}\mathrm{m}}^{-3}}$.

The static stress on the intruder may be modeled using gravitational loading of the grains. In fluids, hydrostatic pressure is given by $\rho g{z}_l$ given gravity $g$ and local plate element depth $z_l$. In granular materials there is a similar depth dependence, however, because the bulk can harbor internal stresses, the behavior differs substantially and ${\sigma }_o=B\rho g{z}_l$ where $B$ can be an order of magnitude greater than in a fluid (Brzinski et al., 2013). In our experiment, the plate extends from the free surface to depth $z$, so $z_l=0.5z$. We assume $B$ is $O(10)$ to estimate that at $z=8$ mm ($z_l=4$ mm) ${\sigma }_o\approx 10\rho g{z}_l=0.06$ Ncm^-2 and at $z=12$ mm ${\sigma }_o\approx 0.09$ Ncm^-2 which are similar to the experimentally measured values.

$B$ can be calculated more precisely using, for example, Coulomb-wedge theory (Guo et al., 2012) or plasticity theory (Kang et al., 2018). However, for the purposes of this study, we were interested in understanding the stress dependence on drag speed and the implications for animal locomotion. Thus, to facilitate comparison with the momentum transfer term, the curves plotted in Figure 4 were calculated using ${\sigma }_o$ estimated from the experimental data. ${\sigma }_o$ was calculated by subtracting $\rho v_d^2$ from ${\sigma }_n$ measured at the three lowest speeds collected at a given $z$ and averaging the result.

Using the empirically obtained static stress, the $\rho v_d^2$ term captured stress as a function of $v_d$ without any free parameters (Figure 4(a)). This was consistent with results for subsurface drag (Brzinski III and Durian, 2010), indicating that, despite the free surface, the physics governing material stresses was similar. While the bulk of the effect of speed was captured, the model was not exact. Future study could explore the origin of this deviation, however, a more thorough investigation of the granular physics is beyond the scope of this paper. Regardless, the predictive power of this simple, time-independent model supported our hypothesis that, across the relevant drag speeds, stress was the result of localized interaction with the material.

Eels in water (a high-Re system) increase the CoM velocity, $v_{CoM}$, by both increasing the speed of wave propagation and altering the waveform to manage the speed-dependent fluid response (Tytell, 2004). C. occipitalis moving on the surface had a linear relationship between the speed of wave propagation and $v_{CoM}$ (Figure 4—figure supplement 1(a)) while the waveform was independent of $v_{CoM}$ (Figure 4—figure supplement 1 (b,c)). This suggested that, despite the $v_d^2$ stress dependence, the aggregate forces responsible for locomotion were not a function of wavespeed but depended only on the pattern of self-deformation. We thus explored the ratio of thrust to drag forces as a function of plate orientation, depth, and speed.

At low-Re, the normal and tangential stresses acting on segments of long, slender swimmers in Newtonian fluids are approximated by ${\sigma }_n=C_n\text{sin}({\beta }_d)$ and ${\sigma }_t=C_t\text{cos}({\beta }_d)$, respectively (Gray, 1955). The constants $C_n$ and $C_t$ are determined by the geometry of a body segment and the viscosity of the surrounding fluid, and the ratio $C_n/C_t$ can be used to approximate swimming performance of a given waveform (Boyle, 2010). Subsurface sand-swimming performance is similarly determined by the balance of thrust to drag, although the prefactors are functions of ${\beta }_d$ (Maladen et al., 2009; Maladen et al., 2011b). Given our hypothesis that the surface slithering snakes were resistive-force-dominated like a low-Re swimmer, we measured the ratio of thrust to drag stress, ${\sigma }_n/{\sigma }_t$. Based on our animal observations and results from the speed-dependent drag experiments, we expected that this ratio depended on drag angle but not drag speed or depth.

We measured stress as a function of ${\beta }_d$ ranging from 0 deg to 90 deg (Figure 3(b)). The ratio ${\sigma }_n/{\sigma }_t$ was largely independent of the drag distance (Figure 5(a)), especially at small ${\beta }_d$. There was a periodic fluctuation of $<1{\sigma }_n/{\sigma }_t$ occurring over several cm appearing in all trials. However, both these fluctuations and the slope of ${\sigma }_n/{\sigma }_t$ as a function of drag distance were small relative to the effect of changing ${\beta }_d$. Thus, we averaged ${\sigma }_n$ and ${\sigma }_t$ to characterize the relationship between plate orientation and stress normal and tangential to its face (Figure 5(b)).

Figure 5

Download asset Open asset

Normal stress increased monotonically with ${\beta }_d$ while tangential stress was approximately constant, gradually falling to zero as ${\beta }_d$ went to 90 deg (Figure 5(b)). The ratio ${\sigma }_n/{\sigma }_t$ as a function of ${\beta }_d$ for both varying speeds and depths collapsed to a characteristic anisotropy curve (Figure 6(a)). The same drag anisotropy curves were measured in poppy seeds and oolite sand and were similarly independent of depth and speed (McInroe et al., 2016). Theappearance of this curve in diverse GM suggests it may be a more general feature of dissipative, deformable materials.

Figure 6

Download asset Open asset

In viscous fluids the constants $C_n$ and $C_t$ are independent of drag angle such that $\frac{{\sigma }_n/{\sigma }_t}{\text{tan}({\beta }_d)}$ is constant. Consistent with results for subsurface drag (Maladen et al., 2009; Maladen et al., 2011b), we find that at the surface $\frac{{\sigma }_n/{\sigma }_t}{\text{tan}({\beta }_d)}$ is a nonlinear function of ${\beta }_d$ (Figure 6(b)). Because tangential stress magnitudes were relatively small and, unlike in viscous fluid, normal stresses rose sharply with ${\beta }_d$ at small angles (Figure 5(b)), thrust and drag forces on the plate were equal (${\sigma }_n/{\sigma }_t=1$) at smaller ${\beta }_d$ than in viscous fluid or subsurface in GM (Figure 6(a)). Reflecting the small angles where the anisotropy was unity, ${\beta }_s$ measured on the snake were small (Figure 6(a), gray histogram).

Plotting anisotropy at constant ${\beta }_d$ as speed and depth varied further illustrated the relatively weak dependence of ${\sigma }_n/{\sigma }_t$ on these variables (Figure 6(c,d)). Especially at small ${\beta }_d$, anisotropy did not depend on $v_d$ (Figure 6(c). ${\sigma }_n/{\sigma }_t$ was more dependent on $z$ at shallow depths, decreasing twofold as $z$ increased from 4 to 8 mm (Figure 6(d)).

The invariance of ${\sigma }_n/{\sigma }_t$ with speed was in accord with our hypothesis of locomotor performance that was dependent solely on the pattern of self-deformation. The efficacy of the simple $\rho v_d^2$ model of the relationship between speed and granular stress indicated that time-dependent structures like vortices were not present, supported by observation of the fast dissipation of disturbance to the grains in both animal and drag experiments. These results indicated that resistive force theory (RFT, Gray, 1955), in particular the same speed-independent granular RFT used to predict subsurface performance in GM (Maladen et al., 2009; Zhang and Goldman, 2014), was a viable candidate for modelling snake motion on the granular surface.

We began by using RFT to estimate the steady-state $v_{CoM}$ of the sand-specialist given the average $v_{seg}=0.67$ ms^-1 measured in experiment, snake length 39 cm, body radius 3.5 mm, and mass of 19 g. We calculated the average $v_{CoM}$ by integrating ${\overrightarrow{v}}_{CoM}$ over a cycle to find the total displacement of the CoM and dividing by the time to complete one cycle (the period). For scale friction equal to that of C. occipitalis, RFT accurately estimated $v_{CoM}$ (Figure 7(a)). Notably, we did not account for material hysteresis in our calculation.

Figure 7 with 1 supplement see all

Download asset Open asset

Subsurface GM swimming performance was improved by reducing scale friction (Sharpe et al., 2015). The surface granular RFT similarly predicted that slithering speed decreased as scale friction increased (Figure 7(a)). This is an intuitive result, as increasing the tangential stress on the side wall by increasing friction of the scales shifts the anisotropy curve down, moving the location of ${\sigma }_n/{\sigma }_t=1$ to larger ${\beta }_d$ (Figure 6(a)) such that force balance requires greater slipping of the body relative to the substrate and thus reduced movement speed.

At the apexes of the body wave, body segments are at ${\beta }_s$ near 0 deg, where ${\sigma }_n/{\sigma }_t<1$ (Figure 6(a)). These portions of the trunk experience drag without contributing thrust. Previous research found snakes moving on firm surfaces lifted these sections of the body off of the substrate (Hu et al., 2009). In the lab we measured 3D kinematics and found that C. occipitalis lifted the wave apexes such that there was a vertical wave on the body which was twice the spatial and temporal frequency of the horizontal wave (Figure 7—figure supplement 1). Therefore, removing these segments from contact with the material reduces drag without a decrease in propulsive force. We included lifting in the RFT calculation by assigning zero force to segments at the wave apexes (Figure 7—figure supplement 1).

Given the low-friction scales of C. occipitalis, including lifting in the calculations did not substantially improve the prediction of animal $v_{CoM}$ (Figure 7(a)). We speculate that these segments may be lifted as a side-effect of the muscle activation responsible for generating the horizontal waveform. It may also be that the lifting is intentional and a buffer against deleterious motor program mistakes or changes in the environment, or the small benefit of lifting these segments is greater than the energetic cost. Because the contribution of both ventral drag and lifting was minor, we did not re-distribute ventral drag forces to account for lifting.

Using the sand-specialist lifting pattern and scale friction we performed test RFT calculations using the drag stresses measured at each depth and found that the effect of depth was less than that of ${\theta }_m$ and less than the experimental uncertainty in $v_{CoM}$ (Figure 7(b)). Further, the peak of the curves did not depend on depth. Therefore, we chose to use force relations measured at a depth of 8 mm. The prediction was similarly insensitive to $v_d$, and captured the relationship between $v_{seg}$ and $v_{CoM}$ measured in the animal (Figure 7(c)).

Despite the observed complexity of motion at the surface–hysteresis, the potential to make and break contacts with the substrate by lifting, and speed-dependence of the drag force magnitudes–the accurate RFT predictions further served to confirm the observation that the animal was in a resistive-force-dominated regime where body inertia was negligible and performance was determined by the pattern and not the speed of self-deformation.

While the total torque about the CoM, summed over all $ds$, at each time was zero, the torque on any given segment was not. We thus used RFT to predict the maximum joint torque, ${\tau }_{m,RFT}$, experienced by any C. occipitalis body segment over an undulation cycle (Figure 8(a)).

Figure 8 with 1 supplement see all

Download asset Open asset

We calculated torque acting on the CoM assuming infinitesimal width and finite segment length. For joint $i$ the internal torque was ${\displaystyle {\tau }_{i,int}=\sum _{k=1}^{i-1}(\overrightarrow{r_k}-\overrightarrow{r_i})\times \overrightarrow{d{F}_k}}$ (Ding et al., 2011). We find this value for each of the segments at all times. ${\tau }_{m,RFT}$ is then the maximum value in this 100 × 70 matrix.

${\tau }_{m,RFT}$ increased as ${\theta }_m$ or $\xi$ decreased. This was both because balancing forces required larger ${\beta }_d$, thus greater force magnitudes (Figure 5(b)), and those shapes ‘stretched out’ the body, creating longer lever arms (see Figure 8(b) and compare shapes in lower left corner to upper right).

The snakes could minimize torque by increasing $\xi$ and/or ${\theta }_m$, however, we did not observe these waveforms (Figure 2). Similarly, snakes did not use waveforms which minimized mechanical cost of transport or distance traveled per cycle (Figure 8—figure supplement 1). We thus hypothesized there were internal factors which ‘penalized’ the low-torque waveforms.

If, for simplicity, we assume the animal moved with no slip regardless of waveform such that the distance traveled per cycle was equal to the wavelength, there were two (not necessarily independent) ways to increase $v_{CoM}$: increase the wavelength by decreasing $\xi$ and/or ${\theta }_m$ (see example waveshapes in Figure 8(b)) or increase the speed of wave propagation (given the linear relationship between propagation speed and $v_{CoM}$, Figure 7(c)).

In limbed organisms, movement speed is related to the interplay between gait parameters like stride length and frequency and physiological concerns like energetic cost (e.g. Hoyt et al., 2006; Kram and Taylor, 1990). We hypothesized that, because the snake performed an escape response in our experiments, their objective was to maximize speed.

For a desired $v_{CoM}$, we explored how quickly a ‘muscle’ segment of fixed nominal length relative to total body length would have to shorten as a function of $\xi$ and ${\theta }_m$. The shortening speed was a function of both how quickly curvature changed along the body for a given waveshape and the frequency needed to achieve the target $v_{CoM}$ for that shape’s stride length. We modeled a body segment as a bending beam with length along the spine $\delta s$ and arclength along the inside of the curve $\delta s^{\prime }$ (see diagram Figure 8(b)). We approximated the amount of shortening the muscles must produce as the difference, $\mathrm{\Delta }s=\delta s-\delta s^{\prime }$, at the point of maximum bending. Using this model we estimated $v_{shorten}=\delta s/\frac{1}{4}T$, the speed the inside of the bend must change length to go from $\delta s^{\prime }=\delta s$ at the inflection points of the wave to the maximum $\delta s^{\prime }$ occurring at the peaks given undulation period $T$ set by the frequency (see Materials and methods for details).

Figure 8(b) is a heatmap of the increase or decrease in the required shortening speed relative to a nominal shortening speed, $v_{short,o}$ (which we take to be that occurring at the mean ${\theta }_m$ and $\xi$ measured in experiment). The actuator speed increased with $\xi$ and ${\theta }_m$. The reason for this was twofold. Firstly, both greater amplitudes and wavenumbers required each segment to enact a greater change in body curvature. Secondly, increasing $\xi$ and/or ${\theta }_m$ decreased stride length, requiring greater frequency to maintain CoM speed. This is somewhat counter intuitive when considering (non-slip) limbless undulatory gaits as compared to limbed gaits. Increasing stride length on a limbed locomotor requires a greater amount of actuation as, e.g., the hip joint must sweep out a larger angle. Therefore to maintain a CoM speed there is a trade off between increasing/decreasing stride length and decreasing/increasing frequency. In contrast, we see here that actuation speed monotonically decreases as the wavelength increases; if we consider only actuation speed, the snakes’ waveforms are under-performing.

This result rationalized why all wave parameters we measured inhabited the center of the wave parameter space (Figure 2); there was a trade off between external forces and internal actuation constraints.

We endeavored to include the trade off between decreasing torque and increasing actuation speed needed to maintain $v_{CoM}$ (now using RFT to estimate the actual stride length given slipping) as ${\theta }_m$ and $\xi$ increased. To do so, we used RFT to calculate the joint-level power of each segment over a cycle for each shape in the ${\theta }_m$, $\xi$ space.

The angle between adjacent segments, $\zeta$, relates to $\theta$ via $\zeta =\frac{d\theta }{ds}$. Joint-level power for each joint at each time was calculated using the rate of change of the joint angle, $P_{i,int}={\overrightarrow{\tau }}_{i,int}\cdot \frac{d{\overrightarrow{\zeta }}_i}{dt}$, where

The power-limited velocity, $v_{pl}$, was the CoM speed for which the peak power generated by any joint over the cycle was equal to a constant peak available power.

In accord with the tradeoff between internal actuation demands and external torques, $v_{pl}$ was maximized in the center of the ${\theta }_m$, $\xi$ space inhabited by the waveform of C. occipitalis (Figure 9(a)). Reflecting the oscillations in $\tau$, this metric had maxima near integer values of $\xi$ (Figure 9(a,b)). Power-limited velocity was maximized at the number of waves (Figure 9(b)) as well as at the attack angle (Figure 9(c)) used by C. occipitalis. We estimated the peak torque output of C. occipitalis’ muscles using dissection of museum specimens (see Materials and methods for details). The attack angle used by the animals was near the point where ${\tau }_{m,RFT}$ increased above the estimated muscle capability. We report both maximal torque output (Figure 9(c) upper horizontal bar) as well as output scaled for the average contraction velocity estimated as the body-wall shortening speed (Figure 9(c) lower horizontal bar). It may be that the individual variation in ${\theta }_m$ reflects differences in peak muscle power capabilities.

Figure 9

Download asset Open asset

Our granular RFT calculations rationalized the stereotyped waveform used by C. occipitalis. The waves of the non-specialists, however, displayed more variation (Figure 2(c)). Using the insights and RFT calculation we developed in our study of the sand-specialist we endeavored to understand how the differences in waveform and morphology impacted performance.

A striking difference between the various species was their aspect ratio, ${\displaystyle A=L/w}$, the ratio of the total length to the diameter of the body at the widest point. A slender snake like Figure 1(b,c) will thus have a higher aspect ratio than a stout one, for example Figure 1(a). We previously discovered that C. occipitalis’ high-aspect-ratio allowed it to move more effectively when swimming subsurface in GM than a low-aspect-ratio sand-swimming lizard (Sharpe et al., 2015).

We characterized the performance of the snakes using the unsigned average slip angle, ${\beta }_s=\text{acos}(|\widehat{t}\cdot \widehat{v}|)$ (Figure 3(a), inset; Sharpe et al., 2015). For each trial we calculated slip along the body at each moment in time and calculated the mean and standard deviation of these values. This variable captures how much the movement of the animal is as if it were moving ‘in a tube’. If the slip is zero we expect each segment follows in the path of the segment before it whereas when slip is large the body is primarily sliding laterally to the midline.

We used RFT to predict the relationship between aspect ratio and slip value. We used the average mass m = 0.76 kg, length L = 0.89 m, ${\theta }_m=48$, and $\xi =2.38$ of all non-sand-specialist snakes with slip < 30 deg and assumed constant body segment length (i.e. all snakes divided into same number of segments), constant body density, $\rho =1000$ kgm^-3, and a circular cross section. The width and length of a snake can then be calculated as ${\displaystyle w=(\frac{m}{\pi \rho A}{)}^{\frac{1}{3}}}$ and ${\displaystyle L=wA}$.

In modeling the various snake species we used frictions 0.1, 0.15, and 0.2. This represents a reasonable range of the friction coefficient between forward-sliding ventral scutes and the substrate during lateral undulation (Hu et al., 2009; Baum et al., 2014), although we note that higher coefficients have been measured on rough material and can be actively modulated by the animal (Marvi and Hu, 2012). We did not account for a difference in intrusion depth which may occur as mass increased. This relationship likely depends on both pressure applied to the GM as well as its properties.

For an animal with fixed mass and internal density, RFT predicted a modest increase in slip with decreasing $L/w$ (Figure 10(a)). Increasing/decreasing scale friction shifted this curve upwards/downwards (Figure 10(a), gray patch). This inverse relationship between aspect ratio and slip was also measured in the animals (Figure 10(a), circle markers).

Figure 10 with 2 supplements see all

Download asset Open asset

There was a transition in the slip versus aspect ratio plot at an aspect ratio of about 26 (Figure 10(a), vertical dashed line). At larger $L/w$ the snakes were, with the exclusion of one exception, moving effectively with an average ${\beta }_s=10.3\pm 3.8$ deg (Figure 10(a), blue crosses). Below $L/w=26$ performance was highly variable. Four of the five species which failed to progress (wave efficiency less than 0.35, Figure 10—figure supplement 1) across the level sand had $L/w<26$ (Figure 10(a), red crosses). The other snakes of $L/w<26$ moved with higher slip on average than their high-aspect ratio counterparts, ${\beta }_s=20.3\pm 5.8$ deg (Figure 10(a), gray crosses).

We found that animals with ${\beta }_s>30$ degrees were those that were ineffective (Figure 10(a), red crosses, Figure 10—figure supplement 1). These animals either did not progress at all (e.g. Figure 1(a,d)), or would have only occasional spurts of forward motion linked by extended periods of undulating in place. Even snakes close to, but less than, 30 deg were able to make consistent forward progress (e.g. Figure 1(e)).

RFT predicted that as attack angle decreased, the slip of a snake with fixed body morphology and wavenumber would increase (Figure 10(b), black line). We examined slip as ${\theta }_m$ changed for those snakes with an aspect ratio less than 26 (and the snake that failed with $L/w$=30). Those which failed to progress (Figure 10(b) red crosses) were at lower attack angles than those which succeeded (Figure 10(b), black crosses)(median and std. of all values ${\theta }_m({\beta }_s\ge 30)=33.2\pm 21.3$ N = 6, n = 4311 measurements, ${\theta }_m({\beta }_s<30,L/w<26)=48.8\pm 24.5$ N = 12, n = 9017, p << 0.01). Those snakes which had $L/w$ of greater than 26 generally used higher attack angles as well (Figure 10 (b, inset) blue crosses).

We used RFT to calculate the performance of individual snakes using that animal’s mass and length and its average wavenumber and attack angle (select points shown in Figure 10(a), square markers). For those snakes with ${\beta }_s<30$ deg and $L/w>26$ the average difference between the measured and predicted values was 5.1 deg (for individuals with more than one trial we compared to the mean of the average slip values, Figure 10(a), black squares). For those with $L/w<26$ and slip $<30$ the calculation was less accurate; the average unsigned difference was 15.3 deg (Figure 10(a), black squares). RFT calculation was the least accurate, average difference 29.9 deg, for those snakes which failed (${\beta }_s>30$ deg, Figure 10(a), red squares). The accuracy of the RFT estimation was inversely correlated with the amount of slip (Figure 10—figure supplement 2).

We observed that the snakes created permanent disturbances in the surface of the GM. Those which moved with low slip created observable piles of sand at the posterior-facing side/s of the body (Figure 1(b,c)), but would not push these piles so far that they would collide with previously-created mounds. In contrast, those which failed appeared to dig themselves into a channel, primarily depositing GM laterally to the long-axis of the body (Figure 1(a)). We noticed that with each undulation these animals would create new piles and then push them until they contacted the mounds of material left behind by previous wave cycles. These material re-interactions occurring during high slip locomotion were not included in the RFT calculations; we hypothesized that once slip was substantial enough that disturbed material was re-encountered in subsequent undulations the material memory impacted performance.

The robot executed serpenoid curves (Equation 1), realized using a Lynxmotion SSC-32 servo controller to dictate the angles between adjacent motors, $\zeta$, (Figure 11(a)). Each motor was offset from its anterior neighbor by a constant phase set by $\xi$. We varied $\xi$ and ${\displaystyle {\theta }_m}$ and recorded kinematics as the robot performed three undulation cycles on the surface of poppy seeds in a 1.2 × 1.8 m² bed. The GM was manually smoothed between trials. Four OptiTrack cameras (Flex 13, Natural Point) tracked infrared reflective markers on the robot (Figure 11(a)). The number of joints on the robot was limited by motor strength. Thus, the robot could only achieve $\xi$ between 1 and 1.4 because lower $\xi$ required more torque than the motors could provide and higher were not well-resolved by the number of joints. The power per unit volume of the motors also limited the aspect ratio; $L/w=11.3$ is lower than that of the stoutest snakes tested.

We characterized performance by measuring slip as in the snakes (Figure 12(a)). We also measured the CoM displacement in a single undulation cycle normalized by the total length (BLC, Figure 12(b)). This variable provided intuition for how effectively the robot was progressing across the substrate.

Figure 12

Download asset Open asset

Performance of the robot was a function of the waveform parameters. This was intuitive given our RFT predictions for the snake, which depended on both $\xi$ and ${\displaystyle {\theta }_m}$ (e.g. Figure 9), and was commensurate with our findings in the biological experiments (Figure 10 (b)). Over the range of $\xi$ available to the robot, performance was more strongly dependent on attack angle than wavenumber, therefore we focused our attention on variations in ${\displaystyle {\theta }_m}$.

We used granular stress relations measured in poppy seeds using a flat plate of the same plastic used to print the robot body (McInroe et al., 2016) in a resistive force theory calculation predicting the robot performance as a function of ${\displaystyle {\theta }_m}$ for fixed $\xi =1$ (Figure 12 (a,b), black curves).

In previous studies, RFT accurately predicted locomotor performance of a number of systems (Zhang and Goldman, 2014) and in our study RFT was accurate for C. occipitalis as well as the successful non-sand-specialist snakes. However, as in the high-slip snakes, RFT under-predicted slip at low attack angles (Figure 12 (a)). Despite the substantial predicted slip, RFT indicated these waveforms would still make forward progress (Figure 12 (b,e)). The robot, however, did not displace at the smallest attack angles tested (Figure 12 (b,c)). Conversely, at high ${\displaystyle {\theta }_m}$ RFT under-predicted displacement of the robot (Figure 12 (b)).

At both low and high attack angles RFT over-predicted the amount of yaw (rotation about the CoM in the horizontal plane) that would occur (Figure 12 (c:f)). The discrepancy was likely because this calculation assumed the robot was always encountering undisturbed material. This assumption worked for subsurface swimming where the GM behaves like a frictional fluid, constantly re-flowing to fill in the spaces cleared by motion of the body. We observed that the motion of the robot deposited some material lateral to the direction of motion. These piles, when re-encountered, appeared to limit yaw. We posited that material hysteresis played an important role in performance.

We observed that material flow was a function of the robot waveform. We recorded overhead video at 120 fps (AOS S-Motion, Baden Daettwil, Switzerland) and estimated grain velocities using particle image velocimetry (PIV, PIVlab ; Thielicke and Stamhuis, 2014a; Thielicke and Thielicke, 2014b). We characterized the overall structure of the flow by measuring ${\psi }_v$, the angle between the velocity vector of the grains at a given point in time and space and the overall direction of motion of the robot (Figure 13(a)). We calculated the probability density to measure a given ${\psi }_v$ over a run for $\xi =1$ and each of the three ${\displaystyle {\theta }_m}$ tested (Figure 13(b)).

Figure 13 with 3 supplements see all

Download asset Open asset

At high attack angles significant amounts of the material were deposited by the posterior-facing body segments, similar to what we observed in the successful snakes (Figure 11(c), Figure 13—video 1). Consistent with this observation, we measured peaks in the probability density of ${\psi }_v$ near ± 180 deg, corresponding to grains which are flowing opposite the average direction of motion of the robot (Figure 13(b), orange curve). When the granular piles moving with the body at later points in time re-enountered these piles the additional stress resisted backward slipping of the robot, increasing displacement above that predicted for a frictional fluid (Figure 12(b)).

The GM being pushed by low ${\theta }_m$ waveforms had a significant velocity component lateral to the direction of motion; as ${\theta }_m$ decreased, the peaks in the ${\psi }_v$ probability density shifted toward ±90 deg (grain flow perpendicular to the average direction of motion, Figure 13(b), green and blue curves). These piles were thus deposited to the side of the robot and with each undulation more material was added to these ‘sidewalls’ which can be seen in Figure 11(b). As the GM being pushed by subsequent undulations encountered these pre-existing walls the robot was apparently unable to overcome the inertia of the grains in the previously created piles, instead depositing more material in the walls without changing their location. The result was that the robot swept out a trough in the GM, completely stalling forward progress and sometimes even moving backward as the trunk ‘rolled’ along the trough walls (Figure 13—video 1). These failures were similar in appearance to those observed in the living snakes (Figure 1(a,d), Figure 13—video 2).

In limbed systems, RFT has often been useful in predicting the performance of the first gait cycle, before material re-interaction could occur (Mazouchova et al., 2013; McInroe et al., 2016). In our limbless robot, however, RFT was incorrect even for the first cycle (Figure 13(c)). Enhancement or degradation of the performance was occurring within the first gait cycle as portions of the body contacted GM disturbed by other segments at a previous time.

This provides some intuition as to why performance of C. occipitalis was robust to variations in the waveform which would have impacted the robot (Figure 2); because this animal’s long body and slick scales yield low-slip slithering, the body and/or piles of GM being pushed by the body are continuously encountering new material through time. The animal is effectively moving in a frictional fluid even though the material is not re-flowing as in subsurface movement.

This raises an interesting subtlety: the tail is always moving through material disturbed by the head, and placing a successful robot waveform back in its own tracks does not decrease performance (Figure 13—video 3). This underlines the idea that locomotive failure does not arise solely because the material is not pristine. Rather, it is the interdependence between the motion of the animal or robot and the evolving material state, both of which depend on the history of the system, that leads to either beneficial or deleterious reaction forces. Despite the complicated interplay between the robot and the substrate, we found the resulting performance was repeatable. For example, for all three attack angles tested at $\xi =1$, each of the three trials for a given waveform pushed material in a similar way and displaced a similar amount in each of the three undulations (Figure 13(b,c)). This indicated that, while the material remodeling was complex, it was deterministic.

Studies found that generalist snakes chose self-deformations in response to available environmental forces (Gray and Lissmann, 1950, Kelly et al., 1997), a strategy which we might expect would be useful on challenging, deformable materials like GM. However, we found evidence that C. occipitalis did not change its waveform when faced with changes to the surroundings (Schiebel et al., 2019) and our robot was executing predetermined waveforms. The ability to sense and respond to induced changes in the substrate is apparently unnecessary given the initial selection of an appropriate pattern of self-deformation.

Although RFT did not consistently predict the robot’s slip, it did provide a heuristic for success and give some intuition as to why a locomotor with the morphology of the robot is challenged by GM. The line in the ${\theta }_m$, $\xi$ space which RFT predicted to be at ${\beta }_s\approx 30$ deg separated the successful and unsuccessful robot waveforms (Figure 14(a)).

Figure 14

Download asset Open asset

We found there were a limited range of waveforms predicted to both have ${\beta }_s\le 30$ deg and be accessible to the robot given its aspect ratio and joint resolution (Figure 14(a)). In contrast, C. occipitalis had a large space of parameters for which slip was small (Figure 14(b)). This suggests that, in contrast to the robot which was sensitive to waveform, the snake was robust to fluctuations in the motor pattern.

The snake also had the benefit of greater resolution than the robot; the animal has ≈ 120 vertebrae (note that in snakes there is variation in vertebral number between species and between individuals of the same species) while the robot has only 10 joints. This flexibility allowed the snake to operate well away from low-performing shapes (low $\xi$ and ${\theta }_m$). Increasing resolution would allow the robot to use more waves on the body, making a greater range of low-slip waveforms available (those in the upper right-hand corner of Figure 14(a)). However, the robot is limited by the size and strength of the motors. Adding more actuators would increase torque on the central joints, preventing accurate realization of the commanded shapes.

• Perrin E Schiebel
• Henry C Astley
• Jennifer M Rieser
• Alex M Hubbard
• Daniel I Goldman

• Jennifer M Rieser
• Perrin E Schiebel
• Henry C Astley

• Perrin E Schiebel

• Perrin E Schiebel
• Henry C Astley
• Jennifer M Rieser
• Christian Hubicki
• Alex M Hubbard
• Kelimar Diaz
• Daniel I Goldman

• Perrin E Schiebel
• Kelimar Diaz
• Daniel I Goldman

• Perrin E Schiebel
• Jennifer M Rieser
• Henry C Astley
• Alex M Hubbard
• Daniel I Goldman

• Perrin E Schiebel
• Jennifer M Rieser
• Henry C Astley
• Daniel I Goldman

• Ken Kamrin
• Shashank Agarwal