Surprising simplicities and syntheses in limbless self-propulsion in sand

In recent years, the study of biomechanics has helped disentangle the roles of control and mechanics during locomotion in natural environments. During animal locomotion, coupled nervous, muscular and mechanical systems produce movements that generate reaction forces to move the body forward. Movement, however, cannot proceed without appropriate environmental interactions. Even common environments like hard ground, water and air display complex interactions with moving animals. For example, analysis of legged locomotion on flat rigid terrain can be complex as a result of impulsive and repeated collisions (e.g. Collins et al., 2005; Holmes et al., 2006). Fluid interactions have been extensively studied (Vogel, 1994; Lauder, 2010), but despite possessing the partial differential equations (PDEs) that describe such interactions (the Navier–Stokes equations), the complexity of fluid flow around deforming bodies makes it difficult to gain a fundamental understanding of locomotion in these environments (Tobalske, 2007; Tytell et al., 2010; Waldrop and Miller, 2015; Gemmell et al., 2016; Cohen et al., 2018; Wise et al., 2018).

In environments composed of granular media (see Glossary), such as sand and loose soil, the situation is seemingly more complicated. Granular media comprise numerous discrete particles which, depending upon loading conditions, can undergo jamming (see Glossary) via internal friction and behave as a solid, or they may slip past each other to yield and flow like a fluid (Gravish et al., 2010). Granular media discussed here will be defined as collections of approximately spherical particles that interact via repulsive contact forces (friction and repulsion). Until recently, we have not had equations for granular materials that describe these fluid- and solid-like interactions under all conditions (see Askari and Kamrin, 2016); however, many insights have been gained through hundreds of years of study (Coulomb, 1776; Schofield and Wroth, 1968; Savage, 1984; Andreotti et al., 2013). For example, provided that stresses do not exceed critical ‘yield’ stresses (see Glossary), granular media can support loads without significant deformation. However, if yield stresses are exceeded, the materials can enter a ‘frictional fluid’ state (see Glossary) that is peculiar compared with true fluids. For example, even at relatively high speeds of movement through granular media, reaction forces are insensitive to rate of movement (Geng and Behringer, 2005; Maladen et al., 2009), in contrast to the strong dependence upon velocity in fluids. This insensitivity is a consequence of the dominance and rate insensitivity of Coulomb friction (see Glossary) between particles and between particle and body elements (Maladen et al., 2009). Further, forces depend on compaction; increasing the packing fraction (see Glossary) by a few percent can lead to a doubling of resistance forces.

Both solid- and fluid-like responses can be important to understand the control of locomotion on granular materials. For example, because the transitions between these two responses are sensitive to the kinematics of the intruding object, animals may modify their movements to promote solidification or fluidization of the medium to improve their locomotor performance (Mazouchova et al., 2010; Marvi et al., 2014). Further, the high frictional forces in sand and the ability of the material to ‘remember’ disturbances (e.g. footprints and tracks) impose additional demands on locomotion. These disturbed regions can have lower yield stress; thus, an animal that interacts with disturbed media may have reduced displacement, setting off a feedback cycle of performance decay that leaves it stuck in its own tracks (Mazouchova et al., 2013; McInroe et al., 2016; Schiebel et al., 2019 preprint).

Despite the challenges associated with the physics of granular environments, many species move on or within these substrates. These include legged (Li et al., 2012; Qian et al., 2015), flippered (Mazouchova et al., 2013; McInroe et al., 2016) and functionally limbless species (Maladen et al., 2009; Sharpe et al., 2015a; Marvi et al., 2014; Astley et al., 2015). In this Review we focus on situations in which experimentation, theory and computation have been integrated to analyze a particular form of locomotion: the use of waves of body undulation to propel the animal on and within dry granular media. Numerous terrestrial vertebrates have elongate, limbless (or functionally limbless) body plans, with most examples from squamate reptiles (Gans, 1975). Functional limblessness has evolved multiple times (Wiens et al., 2006), in deserts and other environments, and is associated with diverse locomotor behaviors (Gray, 1946; Gray and Lissmann, 1950; Jayne, 1986). For example, lateral undulation (‘slithering’) and sidewinding (gaits we will focus on in this Review) both employ posteriorly propagating waves of lateral bending, and sidewinding adds an additional vertical wave of lifting/lowering; lateral undulation and sidewinding are mostly used to move through cluttered substrates or environments dominated by open, granular substrates, respectively (Gray, 1946; Gray and Lissmann, 1950; Jayne, 1986). These behaviors are dominated by frictional forces, with minimal inertial effects (Hu et al., 2009).

Surprisingly, limbless locomotion on granular media provides opportunities to gain insights into these substrates and terrestrial locomotion more broadly; a combination of factors creates a tractable system for examining fundamental concepts in locomotion, such as neuromechanical phase lags, control ‘templates’ (Full and Koditschek, 1999) and morphological and behavioral adaptations. First, dry granular substrates can pose substantial challenges to such locomotors, and thus may drive behavioral or morphological adaptations (Marvi et al., 2014; Tingle et al., 2017). However, unlike more complex terrestrial environments, the granular natural habitat is often homogeneous and obstacles are sparse, making these environments particularly amenable to the creation of laboratory model systems. Tools (Fig. 1A) allowing repeatable and standardized preparation of granular media into ecologically relevant initial conditions are relatively straightforward to develop (in contrast, studies of wet soils and muds often rely on substrate mimics; Dorgan et al., 2005). The rheology (see Glossary) of granular media makes certain locomotion situations relatively straightforward to analyze, and subsurface locomotion is governed by relatively simple forces whose forms are similar to early theoretical approaches taken to understand movement in fluids at low Reynolds numbers (e.g. spermatozoa, nematode worms; Hatton et al., 2013; Goldman, 2014); wet substrate systems display richer but often more complicated interactions (Hosoi and Goldman, 2015; Dorgan, 2015; Sharpe et al., 2015b; McKee et al., 2016). Because environmental forces dominate inertial ones in dry granular systems (as in the world of microscopic organisms), purely kinematic models can capture the motion of animals (Purcell, 1977; Gong et al., 2015; Rieser et al., 2019 preprint). Furthermore, simple mathematical models (‘resistive force theory’) can capture a wide range of granular preparations, from particles of different size, roughness, compaction and incline (Li et al., 2013; Marvi et al., 2014); such models can also apply more broadly to frictional (but non-flowing) interactions (Rieser et al., 2019 preprint).

We now describe the suite of experimental, computational and theoretical tools we have developed to gain insight into granular locomotion, and which can be applied more broadly to other terrestrial locomotion scenarios.

Recent progress in the study of locomotion in granular media has furnished a variety of tools that allow us to investigate and quantify the movements of limbless animals on and within granular media with remarkable thoroughness (Fig. 1). The fluidized bed is one of the most crucial tools – an apparatus that solves the problems of both regional heterogeneities in sand properties and history dependence (e.g. retention of tracks and disturbed media) by allowing the entire sand bed to be reset to a standard state for each experimental trial (Fig. 1A) (Li et al., 2009; Qian et al., 2015). In a fluidized bed, air is propelled by an air blower through a metallic honeycomb supporting a porous plastic surface that is impermeable to granular particles and upon which the medium rests (Fig. 1A). This airflow fluidizes the granular medium, removing any prior disturbances and irregularities. Upon cessation of the airflow, the medium falls into a standardized state (Li et al., 2009). Further modifications to enhance the utility of fluidized beds are available, including shaker motors and air pulses to control the packing fraction of the granular media (Li et al., 2009) and sub-fluidization airflow to reduce yield force (Brzinski and Durian, 2010; Gravish et al., 2010; Zhang et al., 2013). Movement on granular media can be tracked using standard and high-speed videography and motion-capture technology (although observing foot kinematics can be a challenge; Li et al., 2012); subsurface motion can be tracked using cineradiography (Maladen et al., 2009; Ding et al., 2013; Sharpe et al., 2013, 2015a; Fig. 1A).

Finally, physical modeling using robots (i.e. ‘robophysics’) provides a powerful tool to explore the interactions between locomotors and the granular substrate (Aguilar et al., 2016; Fig. 1E). The robophysics approach uses robots to discover principles of locomotion, often through exploration of parameter combinations not observed in animals, thereby determining the consequences of alternative morphologies and kinematic strategies. The repeatability enabled by robots allows more consistent testing and larger sample sizes than with animals, and, when combined with fluidized beds, offers the potential for automated experiments (Qian and Goldman, 2015). Motors with appropriate feedback capabilities can measure force and position data, as well as electrical power consumption, and can provide indirect proxies for measurements that can be difficult to collect from animals (e.g. tendon force-buckles, sonomicrometry, respirometry). Lastly, robots allow validation of the mathematical tools described below (Maladen et al., 2011a,b).

Small systems of granular media can be modeled by computers via ‘discrete-element modeling’ (Ding et al., 2012; Maladen et al., 2011a,b; Fig. 1B), in which every particle's interactions with other particles and with intruders (i.e. objects moving through the medium) are simulated. Although this approach is accurate, it is also computationally intensive even for systems of millions of particles; real sand beds can contain many orders of magnitude more grains. The ability to understand locomotion in the context of granular media has been substantially increased via the mathematical tool of resistive force theory (RFT) (Box 1; Maladen et al., 2009, 2011a,b; Li et al., 2013), originally developed to estimate the forces on low-Reynolds number (Re) swimmers in Stokes' flow (Gray and Hancock, 1955; Fig. 1C,D). RFT assumes that the total force acting on an object, such as an animal's body, is a linear, independent summation of the forces acting on its constituent parts (Box 1). The granular drag forces acting on body segments are measured in systematic experiments using intruding objects of various sizes, shapes and orientations moving though granular media in fluidized beds (Box 1). RFT calculations using these empirical measurements capture crucial aspects of movement through dry granular media, especially if history-dependent effects (e.g. footprints) are minimal and intruders move at moderate to slow speeds (Box 1) (Marvi et al., 2011; Mazouchova et al., 2013; Aguilar et al., 2016; McInroe et al., 2016).

Recent work has revealed that granular RFT emerges from features of a set of PDEs known as ‘frictional plasticity models’ (Askari and Kamrin, 2016; S. Agarwal, A. Karsai, D. I. Goldman and K. Kamrin, in preparation). Until recently, we have not had such equations; thus, this can be considered a major advance with great potential for future analyses. Solving the plasticity PDEs is computationally intensive; however, because RFT is accurate and straightforward, we can apply this method to calculate forces for infinitesimally small segments of a moving body (Li et al., 2013), and subsequently predict the performance of hypothetical waveforms (Maladen et al., 2009, 2011a,b) as well as power (Ding et al., 2012) and joint torques (Ding et al., 2012; Sharpe et al., 2013; Fig. 1C,D). Thus, as in the robophysical studies, the use of RFT in modeling allows the exploration of movement patterns beyond those observed in animals. Likewise, we can test whether the animals are indeed moving optimally, and determine the consequences of alternative movement patterns (Fig. 1C,D). Further, combining RFT with analytic tools from geometric mechanics (Hatton et al., 2013) enables visualization of optimal or alternative gaits and limb-coordination patterns, subject to kinematic constraints and body morphology, without the need for detailed calculations (Fig. 1F).

Our comparative approach uses the mechanical properties of sand, the diversity of morphology and behavior of squamate reptiles, and mathematical tools such as RFT and geometric mechanics to describe and test the functional principles of locomotion. This robust approach can test basic principles of locomotor performance, with clear applications to evolutionary biology, biomechanics and, ultimately, robotics. Furthermore, our approach allows bi-directional transfer of knowledge between biology and physics/robotics: testing of observed and quantified biomechanical performance by animals can lead to physical insights, which in turn can allow exploration of parameter ranges beyond those observed in biology, thereby allowing direct comparison between observed behaviors and hypothetical alternatives. This reciprocal approach informs robotic potential and allows for systematic testing of biological systems, which facilitates objective tests of adaptation and optimization, and even allows us to understand the influences of phylogeny with respect to locomotor opportunity or constraint.

In order to demonstrate the potential insights to be gained from the study of locomotion in granular media, below we discuss recent work focused on three species that specialize in movement on or through sand: the sandfish skink (Scincus scincus), the Mohave shovelnose snake (Chionactis occipitalis) and the sidewinder rattlesnake (Crotalus cerastes) (Fig. 2A–C). These animals are ‘sand specialists’, because they not only inhabit dune fields but also have developed specialized modes of locomotion: sand swimming in S. scincus and Ch. occipitalis, and sidewinding in Cr. cerastes. These case studies will illustrate how application of the tools and methods for granular systems can provide insights into physiology, behavior, control and evolution.

The sandfish (S. scincus) (Fig. 2A) is a quadrupedal skink, native to North Africa and the Arabian Peninsula; it takes its name from a remarkable escape behavior of rapidly diving head-first into the sand and then submerging (Arnold, 1995). Although the sandfish possesses prominent limbs, these adduct completely during sub-surface burial, resulting in functionally limbless undulatory propulsion (Maladen et al., 2009). Cineradiography of subsurface locomotion in a fluidized bed reveals that sandfish achieve rapid subsurface locomotion (∼1.2 body lengths s⁻¹, up to 4 Hz), but the body segments show significant rearward slipping due to media yielding around the body (Maladen et al., 2009). To facilitate analysis, body posture was characterized using the number of waves on the body (ξ) and relative body curvature (κ_mλ_s, the arc length of one wave divided by the maximum radius of curvature; Sharpe et al., 2015b), a non-dimensional metric for the degree of bending, which is independent of both animal size and number of body waves. Performance is described using metrics such as estimated center-of-mass velocity, wave efficiency and slip angle (Maladen et al., 2009; Sharpe et al., 2015a). Wave efficiency (η – forward body speed divided by wave propagation speed) or undulation efficiency (η_u, forward distance moved per cycle divided by the arc length of one wave) quantify how effectively body motion is converted to forwards motion, with no slip resulting in a value of η=1 and no forward motion resulting in a value of η=0 (Maladen et al., 2009; Sharpe et al., 2015b). Slip angle (β_s) is defined as the angle between the vector tangent to the body axis at a given point and the velocity vector of that same point (Fig. 1C); if there is no slip, each point follows the path of the points before it and β_s will be zero, whereas any slipping will produce higher values. Unlike η, β_s can be calculated within a cycle and along the body, allowing analysis of where the slip occurs along the body and in which phase of the cycle, though it can also be averaged to generate a whole-animal measurement (Sharpe et al., 2015a).

Slip in sand swimming causes a decrement to locomotor performance, raising the question of whether such slip is unavoidable, and, if so, whether the sandfish are using the waveform with the maximum undulatory efficiency (η_u) possible for them. Sandfish use only a limited range of waveforms during sand swimming (Maladen et al., 2009), so exploration of alternative kinematic strategies requires a combination of mathematical and robotic modeling. Applying RFT to simulated sandfish waveforms allows calculation of η_u across a range of body configurations, including some not observed in the animals, showing that the sandfish operates near the optimum waveform for its morphology in terms of η (Maladen et al., 2009), η_u, β_s (Sharpe et al., 2015a) and mechanical cost of transport (mCoT; Sharpe et al., 2015a; Fig. 3A–C). These calculations were tested using a robotic model, and the results closely match the robot-specific RFT results (Maladen et al., 2011a,b), as well as results of discrete element modeling (Ding et al., 2012; Maladen et al., 2011a,b). Results from the robotic model also support the fundamental relationship seen in the calculations involving sandfish. Torque calculated using RFT for different segments within the body closely matches observed electromyographic signals from burrowing sandfish (Sharpe et al., 2013; Fig. 3D). This provides a potential explanation for the generalized phenomenon of neuromechanical phase lag (the temporal delay observed between muscle activity and the resultant axial bending in undulatory aquatic locomotion across many species) (Ding et al., 2013).

The efficacy and efficiency of limbless locomotion in granular systems, both on the surface and beneath, raises the question of why sandfish retain their prominent legs. During the burial phase of sand swimming, sandfish fold their legs back against the body much like a swimming crocodile (Maladen et al., 2009). Yet, limb reduction and loss is common among scincid lizards (Wiens et al., 2006; Whiting et al., 2003), highlighting an apparently exceptional limbed morphology in the sandfish. To examine this, the motion of the sandfish limbs was tracked prior to, during and after burial, and the effect of limb restraint on burial speed was investigated (Sharpe, 2013). Sandfish use all four limbs to propel themselves on the surface and to plunge their snout into the substrate, but each limb pair folds back against the body as soon as it is submerged (Fig. S1A–C) (Sharpe, 2013). Temporarily restraining the limbs in various combinations revealed the crucial role of the limbs (particularly the forelimbs) in sand burial (Fig. S1D–F) (Sharpe, 2013). Thus, the limbs are used on the surface, and appear to be crucial for crossing the surface–substrate boundary.

The Mohave shovelnose snake (Ch. occipitalis) independently evolved sand-swimming behavior, thereby allowing for a comparison with the sandfish to test the generality of these principles. Its drastically different morphology allows further exploration into the mechanics of subsurface locomotion in sand (Sharpe et al., 2015b). Chionactis occipitalis outperforms sandfish on every metric (β_s, η and mCoT) except speed, with burrowing speed 1/10th that of the sandfish (Sharpe et al., 2015a) (Fig. 3A–C). RFT modeling shows that the improvement in snake performance is due to both lower friction between the scales and the sand and the more elongate, flexible body form, which enables the snake to access combinations of relative body curvature and number of body waves that are anatomically inaccessible to the sandfish (Sharpe et al., 2015a; Fig. 3A–C). However, both species perform close to the optimal kinematics for their body form (Sharpe et al., 2015a; Fig. 3A–C), potentially due to strong evolutionary selective pressure to minimize the tremendous mechanical demands of sand swimming. This is an insight only possible through the combined use of living animals, cineradiography, fluidized beds, RFT and robotics in a granular system.

Limbless vertebrates moving on the surface of granular media use all known methods of limbless vertebrate locomotion (rectilinear, concertina, lateral undulation and sidewinding; Table S1), but here we focus on sidewinding. Sidewinding illustrates a crucial aspect of movement on the surface of granular media: the ability to control the forces (and minimize yielding) on body segments by lifting and lowering portions of the body (Jayne, 1986; Marvi et al., 2014).

Sidewinding involves undulation in both the lateral plane (typical of snakes) and the vertical plane (Fig. 2C, Fig. 4A,B). It is known in a few species of viperid snakes moving across dry sand (Mosauer, 1932; Brain, 1960; Jayne, 1986; Gans and Kim, 1992), and is also used by homalopsid and natricine snake species to traverse wet granular media (Jayne, 1986). Unlike concertina and lateral undulation, which make sliding contact with the ground that may cause yielding of granular substrates, sidewinding snakes maintain approximately (low slip) static contact points with the sand (and rigid substrates), lifting and lowering the body segments as they move (Fig. 2C, Fig. 4A,B; Rieser et al., 2019 preprint).

The net motion of the snake that emerges from these lifting and lowering patterns resembles that of a flattened ‘helical spring’ (Mosauer, 1930) or ‘virtual tread’ (Hatton and Choset, 2010; Gong et al., 2012) rolling across the ground, although the snake's dorsal surface remains upwards at all points unlike the spring. This motion can be decomposed into two orthogonal body waves, vertical and horizontal, ±π/2 (±90 deg) out of phase with each other, resulting in posteriorly propagating waves of lowered static contact and lifted movement, such that a point along the body is cyclically lifted clear of the sand, moved forward, then placed back into static contact (Fig. 2C, Fig. 4A,B) (Marvi et al., 2014; Astley et al., 2015). This two-wave template allows not only easier understanding of a visually confusing locomotor mode (Pope, 1955) but also the modeling of sidewinding in snake robots with vertical and horizontal degrees of freedom, either by application of the neuromechanical template (see Glossary; Marvi et al., 2014; Astley et al., 2015) or by two-wave-based mathematical analyses of motion-capture data to gain insight into template modulation to produce complex behaviors (Gong et al., 2015) (Fig. 4B).

Control of contacts with the substrate is essential to sidewinding, both in generating the fundamental locomotor template and in the modifications of this template to allow effective locomotion under different conditions. Lifting allows some body segments to interact with the substrate while others do not, a trait also used in lateral undulation (Hu et al., 2009), albeit less dramatically. In the absence of the vertical wave in the template, sidewinding simply is lateral undulation (Astley et al., 2015) – the plesiomorphic locomotor behavior in squamates, legged or non-legged. This vertical component of the template allows for subtle control, such as the ability to control the relative length of the body segments that are in contact with the granular media, which is crucial in order to ascend inclined sand with minimal slip in spite of the decreasing granular yield forces of the sand with increased incline (Fig. 4D; Marvi et al., 2014). This strategy has been replicated in a robotic model to enable successful climbing of granular slopes, and the robophysics approach has been used to test combinations of contact length and incline not observed in the snakes, thus determining the consequences of parameters not observed in nature (Fig. 4F) (Marvi et al., 2014). These robophysical tests revealed that snakes are near the optimum of a range of effective kinematics; this range decreases with increasing incline (Fig. 4F; Marvi et al., 2014).

Sidewinding allows snakes to overcome the challenges of surface locomotion on granular media over a wide range of conditions (Marvi et al., 2014; Astley et al., 2015); however, clear, cyclic sidewinding is known only from a handful of caenophidian species (Brain, 1960; Gans and Mendelssohn, 1972; Jayne, 1986), mostly vipers (Gans and Mendelssohn, 1972), with the additional description of a similar, but non-cyclical behavior in a green anaconda (Eunectes murinus; Ryerson and Horwitz, 2014). The taxonomic distribution of these species across multiple clades of alethinophidian snakes raises the possibility that either sidewinding has convergently evolved multiple times or, more parsimoniously, that virtually all snakes can sidewind, to some extent, if presented with suitably challenging granular substrates.

To test these alternatives, we encouraged 30 snake species from a range of taxonomic groups and habitats in the Zoo Atlanta collection (Table S1) to move on our fluidized bed at 0, 10, 20 and 27 deg inclines, the last being the avalanche angle of the sand. Only two of 30 species used sidewinding in any condition (Cr. cerastes and the Northern watersnake, Nerodia sipedon, the latter only very briefly), whereas the remaining species never performed sidewinding even when failing to move using other modes (Table S1, Movie 1; Marvi et al., 2014). Instead, most species used either lateral undulation or concertina locomotion, sometimes resorting to the latter when the former failed, and none were capable of moving uphill at the avalanche angle like the sidewinder. Although strictly proving the absence of a behavioral trait is impossible, it is unlikely that failing animals would refuse to use an effective alternative locomotor mode if it was available to them. This suggests that sidewinding is not present in most snakes and therefore evolved independently several times.

The strongest predictor of performance on sand was taxonomy, specifically whether the snake was a viper – 25% of vipers we tested failed to move even on level sand, whereas all non-vipers succeeded in doing so, and only two vipers could ascend granular inclines (Cr. cerastes and the speckled rattlesnake, Crotalus mitchellii, the latter relying on a mix of rectilinear and concertina locomotion; Table S1). Although the mechanistic explanation for this discrepancy remains unknown (and is the subject of ongoing investigation), it suggests that the prevalence of sidewinding in vipers is not due to some inherent advantage of vipers (Gans and Mendelssohn, 1972), but rather is due to a fundamental locomotor deficit that prevents vipers from accessing sandy habitats without the evolution of sidewinding. This matches observations of viperid species distributions across the Mojave Desert, where several viper species inhabit consolidated terra firma substrates or rocky slopes, but regional sand dune systems are exclusively populated by sidewinders (J.R.M., personal observation).

Among the snakes using lateral undulation on the granular media, some individuals did not make any progress (e.g. Fig. 5A,D), some experienced high slip but made effective forward progress (Fig. 5E,G), whereas others, including the desert specialist Ch. occipitalis, moved with little slipping (Fig. 5B,C,F,H). Effectiveness was related to how the animal's motion remodeled the granular substrate – the tracks. Snakes that failed pushed material lateral to the long axis of the body (Fig. 5A) whereas successful snakes pushed material both laterally and posteriorly (Fig. 5B,C).

Modeling lateral undulation on the surface is complicated by the material hysteresis. However, RFT accurately captured the performance of snakes that did not slip enough to re-encounter their own tracks, and it predicted trends, giving insight into anatomic and kinematic features impacting performance. As in subsurface sand swimming, RFT predicted that performance depended on body elongation (Fig. 5I, black curve; L/w is total body length divided by the width at the widest part; Fig. 5I, inset) and kinematics (Fig. 5J, black curve). In accordance with this prediction, long, slender snakes moved with low slip (Fig. 5I, blue and green markers), with one exception (a viper, Bothriechis schlegelii), while slip was high among stout snakes. Some stout snakes were still able to progress (e.g. Fig. 5E), although slip was variable (Fig. 5I, gray markers), whereas others failed to move, undulating in place as in Fig. 5A (Fig. 5I, red markers). Among stout snakes, those using high relative curvature waveforms moved with lower slip than those at lower relative curvature (Fig. 5J), a trend captured by RFT (Fig. 5J, black curve). RFT predicted performance of those snakes which succeeded (Fig. 5J, black squares), but underpredicted slip of snakes which failed (Fig. 5J, red squares). This was because the model did not include material remodeling; high-slip snakes re-encountered previously disturbed material, thus experiencing forces dependent on the changed substrate state.

Although RFT can calculate the consequences of a particular waveform, exploration of alternative behaviors is time consuming and non-intuitive. However, by combining RFT with geometric mechanics, we can understand the performance consequences of a wide range of parameters in an intuitive way, fostering insights into the mechanical and evolutionary consequences and tradeoffs of behavior. The minimal role of inertia in many forms of movement on or in dry granular media has facilitated the application of this mathematical framework for locomotion (referred to as ‘self-propulsion’) developed decades ago in physics. This framework can unify our understanding of locomotion by allowing diagrammatic visualization of effective (and ineffective) patterns of body undulation and limb motion. The motion of microorganisms in fluids at low Reynolds number (Purcell, 1977) inspired the creation of a systematic approach to link ‘self-deformations’ (changes in body shape and limb position) to translations and rotations in the environment (Wilczek and Shapere, 1989; Murray and Sastry, 1993; Walsh and Sastry, 1995; Ostrowski and Burdick, 1998; Melli et al., 2006; Shammas et al., 2007; Morgansen et al., 2007; Avron and Raz, 2008; Hatton and Choset, 2011, 2013, 2015). Although this geometric framework was originally developed for low-Re swimmers, robotic models showed it could be applied in granular media (Hatton et al., 2013), where the highly damped ‘frictional fluid’ aspects of the granular interaction result in motion that depends largely on sequences of shape changes (thus the term ‘geometric’ mechanics).

The essential approach of geometric mechanics is to analyze the effects of transit around closed (cyclic) paths in a ‘shape space’. Such a space contains representations of all possible shapes (postures) of the organism (see examples in Fig. 6A–C), and short segments in a path in the space connect different instantaneous shapes via ‘self-deformations’. Broadly speaking, if we determine the coupling between small self-deformations and real world displacements, we can calculate the net movement over a cycle by summing the small displacements. This sounds simple, but solving the inverse problem of which cycles best combine the available displacements is computationally expensive and can be difficult to visualize. In recent years, advances have solved some of these difficulties (Hatton and Choset, 2015). Although initially limited to simple shapes with few joints, recent developments have allowed generalization to undulating, continuous swimmers in granular media (Hatton and Choset, 2011; Gong et al., 2016; Ramasamy and Hatton, 2016, 2017; Rieser et al., 2019 preprint), by representing these body shapes as a combination of two waves whose curvature varies sinusoidally in position along the body (called serpenoid waves; Umetani and Hirose, 1975); this system can be generalized to more than two shape variables if needed (Ramasamy and Hatton, 2019; Bittner et al., 2018). A geometric mechanics MATLAB software package is available at https://github.com/OSU-LRAM/GeometricSystemPlotter/releases (Ramasamy and Hatton, 2019).

For systems with two shape parameters (e.g. two joint angles, or the amplitude of two modes of continuous bending), the results from the geometric mechanics approach can be visualized as a set of scalar ‘constraint curvature functions’ (which we will refer to as CCF; see Glossary) over the shape parameter space (with one CCF for each direction of translation or rotation); Fig. 6A,B shows forward displacement CCFs for the sandfish and shovelnose snake. The beauty of the approach is that, given a path in the CCF space (again, a pattern of cyclic self-deformation), we can immediately determine whether this is a ‘good’ or ‘bad’ gait (in terms of displacement) simply by summing the amount of positive and negative surface enclosed (red and black colors in Fig. 6A,B). Given that the CCFs are generated from the environmental interactions and body morphologies, these tools therefore allow evaluation of the observed behavior of an animal and determine whether the animal is optimizing for a given performance metric (Nishikawa et al., 2007; Marvi et al., 2014; Astley et al., 2015; McInroe et al., 2016), as well as predicting the consequences of alternative morphologies and behaviors that are not observed in the animals.

Geometric mechanics analysis can also be used for locomotion on the surface of the sand, including locomotion with discontinuous contacts, as in the case of sidewinder rattlesnakes. Mapping the contact region onto the posture shape space (Fig. 6D) allows calculation of a CCF for sidewinding (Fig. 6E). These calculations can be repeated for animals moving on rigid surfaces with both high and low friction and compared with previously collected kinematic data (Schiebel et al., 2019 preprint); this shows that animals use 1.5 waves per body length. Track angle can be calculated from body waveform, allowing comparison with many more species (Rieser et al., 2019 preprint) (Fig. 6F, inset). Observed data and geometric mechanics calculations agree well, showing that snakes perform better on rigid surfaces (possibly analogous to hard-packed sand in portions of their habitat), but perform close to optimally under both conditions (Fig. 6).

Thus, we can return to the performance data for the species discussed above (e.g. Fig. 3) with a new ‘geometric’ eye on the problem. If we make the approximation that the animals use traveling body waves of sinusoidal local curvature, this leads to circular paths in the CCF composed of two modes that approximate the animals' shape changes (though the actual path of the animal may be non-circular). Remarkably, the animals operate near paths that (out of all circular paths) achieve most displacement per undulation cycle (Fig. 6C), those that essentially follow zero-contours (white regions in Fig. 6A,B); thus, we can conclude that these traveling wave patterns enable the observed rapid movement within and on granular material.

The great power of the geometric diagrams is to provide insights into the functional consequences of an animal's choice of a particular path through space, potentially allowing us to understand why an animal does or does not choose a particular path. That is, geometric mechanics provides candidate control templates and facilitates our understanding of how these can/should change as a function of environment, morphology, body properties, etc. For example, because muscular contractions have a limited peak power and because speed is the product of the displacement per cycle and the cycle rate, further study is needed to determine how muscle physiology, morphology and interactions with the surrounding media influence the effort required for different cycle geometries. The results of these studies will allow us to geometrically identify speed-optimal gaits that balance displacement per cycle against cost per cycle (Ramasamy and Hatton, 2016). Interestingly, these tools also work for the sand swimmers and the sidewinders, indicating their applicability to both fluid-like and solid-like situations. Future research directions should thus include overlaying models of muscle dynamics onto the geometric locomotion framework, as well as collecting new in vivo data to test whether snakes may in fact be using a non-uniform wave amplitude (which would be an oblong path in Fig. 6B) to better take advantage of the geometry of the CCF.

These examples have shown how the study of locomotion on granular media can give substantial insights into the fundamentals of limbless locomotion, including such concepts as neuromechanical phase lags, control templates and morphological and behavioral evolution (Ding et al., 2013; Sharpe et al., 2013, 2015a; Marvi et al., 2014; Astley et al., 2015). Counterintuitively, the homogeneous and speed-independent, yet mechanically complex and demanding nature of these substrates makes mathematical analysis and robotic modeling more feasible than for systems such as fluids, while still imposing sufficient demands on locomotion to prompt interesting and ecologically relevant tests of evolutionary adaptations in nature. Furthermore, the triad of biology, physics and robotics applied in a coordinated manner allows the exploration of alternative morphologies and behaviors that are not testable in nature, in order to determine whether the animals are truly optimizing performance and the consequences of alternative character states. Each of these tools has strengths and weaknesses, but when used in combination, each can overcome the limitations of the other, allowing for insights that are both deeper and more broadly applicable across species. This approach thus enables the construction of hypothetical mechanical analogs to fitness landscapes, quantifying the performance consequences of both observed behaviors and hypothetical alternatives in order to visualize tradeoffs and both local and global optima. Considering the advances enabled by this method, we suggest that granular systems and sand-dwelling animals can function as a ‘model system’, much like model taxa, to enable deeper understanding of locomotor behavior and control in challenging environments.

To fully realize the potential of such systems, we must address several challenges in this area. First and foremost are the consequences of disturbance of the media by the animal itself, leaving tracks with which it may interact during subsequent locomotor cycles (Mazouchova et al., 2013; Schiebel et al., 2019preprint), or disturbances from other animals or forces in nature. Although such disturbances are readily observable and are certainly the norm in nature, our ability to model these disturbances is in its infancy (Mazouchova et al., 2013; Schiebel et al., 2019 preprint). Recent work indicates the existence of a predictive model for the impact of material remodeling on the granular reaction forces and thus locomotor performance (Schiebel et al., 2019 preprint). Such a model could be incorporated into a granular-state-dependent RFT that would elucidate how organisms can mitigate deleterious effects of material disturbances or even take advantage of remodeling to further improve performance (Schiebel et al., 2019 preprint).

Second, the robotic and mathematical models discussed are ‘open loop’ in design, without sensory feedback or behavioral modulation based on the environment. Although this simplification is useful for understanding robust behaviors that maintain performance without feedback, a deeper understanding of animal locomotion must incorporate the feedback and control mechanisms present in real animals. Our approach has been useful to discover control templates (Full and Koditschek, 1999), and we posit that the granular media systems will facilitate the discovery of anchors – models that describe how animal physiology and morphology integrate to produce the template behavior. As this field develops, we look forward to solutions to both of these shortcomings, with the commensurate increased understanding of these fascinating systems. Principles discovered in this tractable, yet complex locomotor substrate may provide insights into more complex substrates, such as heterogeneous and cohesive granular media (Winter et al., 2014; Qian and Goldman, 2015; Dorgan, 2015).

The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the US Government. The US Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.