Limbless animals like snakes inhabit most terrestrial environments, generating thrust to overcome drag on the elongate body via contacts with heterogeneities. The complex body postures of some snakes and the unknown physics of most terrestrial materials frustrates understanding of strategies for effective locomotion. As a result, little is known about how limbless animals contend with unplanned obstacle contacts. We studied a desert snake, Chionactis occipitalis, which uses a stereotyped head-totail traveling wave to move quickly on homogeneous sand. In laboratory experiments, we challenged snakes to move across a uniform substrate and through a regular array of force-sensitive posts. The snakes were reoriented by the array in a manner reminiscent of the matter-wave diffraction of subatomic particles. Force patterns indicated the animals did not change their selfdeformation pattern to avoid or grab the posts. A model using open-loop control incorporating previously described snake muscle activation patterns and body-buckling dynamics reproduced the observed patterns, suggesting a similar control strategy may be used by the animals. Our results reveal how passive dynamics can benefit limbless locomotors by allowing robust transit in heterogeneous environments with minimal sensing. snake | locomotion | neuromechanics | biomechanics | collisions D isentangling the roles of environmental mechanics, physi- ology, and neurology to discover principles of movement in natural environments requires integrating insights from disciplines including neurobiology, biomechanics, control theory, and soft-matter physics (1)(2)(3)(4). While progress in such integration has been made in locomotion in homogeneous environments [open fluids, flat hard ground, and dry sand (4-6)], understanding movement in heterogeneous terrain remains a frontier in locomotion studies (7)(8)(9). This is in part because interactions in these environments can change discontinuously and unexpectedly, making it unclear how biological control schemes and bodyplans accommodate unplanned collisions.
In organismal neuromechanics (10), locomotion control is classified on a spectrum between a closed and open loop (11), each dealing with heterogeneities in different ways. We will consider these defined as the relationship between information about the surrounding terrain and consequent determination of the self-deformation pattern. During the former (typically associated with careful, deliberate movements) the animal uses sensory input to collect information about the terrain and self-deforms in response to novel interactions (12,13). In contrast, rapidly moving animals (9,14) can rely on the passive dynamics of mechanical structures to rapidly reject unexpected perturbations without additional input from the nervous system.
While the strategies used by limbed animals to contend with collisions have been carefully studied (9,14), little is known about how body-undulating organisms like nematode worms, eels, and snakes handle unplanned interactions (e.g., turbulence or obstacles) (7,15,16). Understanding movement strategies in heterogeneous environments is particularly difficult in terrestrial animals like snakes, which often rely on complex terrain heterogeneities to propel themselves. The only gait shared by all limbless, elon-gate vertebrates is lateral undulation (17). This is the commonly seen slithering motion in which primarily planar curves of the trunk press against heterogeneities-whether discrete obstacles (16), those created and subsequently used by the body bends such as piles of sand (18), or frictional anisotropy introduced by the structure of the integument (19)-to generate propulsion.
Coordinating interaction between the many degree-of-freedom trunk and complex terrain is a nontrivial task as evidenced by the failure of snakes challenged to move in a novel terrain (20). It is unclear how limbless organisms contend with unexpected collisions and whether passive mechanics could facilitate robust transit as in limbed systems. Previous studies of generalist snakes (those with a diverse geographical range encompassing a variety of habitats) found that the waveform was related to the density of heterogeneities (21), and subtle, local deformation was precipitated to use obstacles (16). These complex shapes confound understanding of the relationship between terrain, shape change, and performance.
Here, we take a first step to understanding the neuromechanics of laterally undulating snakes in heterogeneous environments
Snakes inhabit environments composed of heterogeneous materials, controlling their body-terrain interactions to generate propulsion. Such complexity makes it challenging to understand the interplay of body mechanics and neural control during obstacle collisions. To simplify, we studied a desertdwelling snake with a stereotyped waveform moving in a laboratory heterogeneous terrain, an array of posts embedded in a sand-mimic substrate. Compilation of hundreds of trials revealed multipeaked "scattering" patterns, reminiscent of diffraction of subatomic particles. A model incorporating muscle activation patterns and body buckling recovered the mechanical diffraction pattern, indicating passive dynamics facilitates obstacle negotiation without additional neural input. Our study demonstrates the importance of mechanics in snake locomotion as well as the rich dynamics in collisions of self-propelled systems.
PHYSIOLOGY by working with a relatively simple system, a desert-dwelling snake that relies on a stereotyped self-deformation pattern to move within its habitat, composed largely of homogeneous sand but containing sparse obstacles (SI Appendix, Fig. S1). In the absence of the ability to interrogate the motor control system in freely moving snakes (unlike the increasing number of tools in microscopic swimmers, e.g., ref. 22), we explored the animals' response to unexpected terrain interactions using a "scattering" approach (23)-studying the kinematic and dynamic outcomes of collisions with heterogeneities. Determining principles of body coordination in complex terrains will help improve mobility of snake-like robots (20,24).
Biological Model. The shovel-nosed snake, Chionactis occipitalis (Fig. 1A), transits open desert between the cover of larger flora to forage or escape threats (18). This intervening terrain consists The first two components, PC 1 and PC 2 , captured 86% of the variance and were well-fit by sinusoids (dashed gray lines; PC 1 R 2 = 0.98 and PC 2 R 2 = 0.99) approximately π 2 out of phase (0.45π ± 0.04π). (F) κ can be approximated by κ(s, t) ≈ α 1 (t)PC 1 (s) + α 2 (t)PC 2 (s). A plot of α 1 (t) versus α 2 (t), colored by frame number, revealed trajectories moved clockwise along a circular path which, combined with the PCs, produces a traveling sinusoidal wave.
of a sand substrate interspersed with sparse heterogeneities such as small plants and twigs (Fig. 1B). This desert-dwelling species uses a stereotyped traveling wave when moving quickly (30 to 80 centimeters per second, 0.8 to 2.1 body lengths per second) on the surface of sand (18) (Movie S1).
Fig. 1C is an archetypal example of the sinuous waveform used by this species when moving on spatially uniform and yielding substrates like sand or the sand-mimic substrate, high-pile carpet, used in this study to increase the rate of data collection (SI Appendix). This traveling wave initiates near the head and passes posteriorly with little variation along the arclength, s, of either maximum curvature, κm, or wavelength (Fig. 1D). Principal component analysis (PCA) (inspired by ref. 25) is a method of dimensionality reduction in which one calculates the eigenvectors of the curvature covariance matrix. The resulting principal components (PCs) form an orthogonal basis describing the variation of κ along the body. PCA revealed that κ was well-approximated by a serpenoid curve (26) κ(s, t) ≈ κmsin( 2πL ξs + ωt) with κm = 25.2 ± 3.0 m -1 and wavenumber ξ = 2.0 ± 0.3 (Fig. 1 E and F) (27) (the waveform changes in time according to ωt; however, because in our system inertia is dominated by damping, we will not include ω in analysis and ω = 1 in all models).
We previously discovered that targeting this serpenoid waveform conferred benefits to the snake's movement on the surface of homogeneous sand (28), providing rationale for the conserved appearance of the waveform across individuals and trials. Given these locomotor benefits we next examined how the wave pattern changed upon unexpected collisions. Would the animal actively alter the waveform to either avoid or use reaction forces from the posts or could it use passive dynamics to transit the array without active changes to the serpenoid wave?
We modeled the sand and sparse heterogeneities of the desert terrain as a row of six rigid, force-sensitive posts (Materials and Methods) embedded in the carpet substrate (Fig. 2 A and B). We challenged C. occipitalis (N = 8) to travel across the substrate and through the post array. High-speed video captured kinematics and custom MATLAB software digitized the snake midlines for analysis (27,29) (SI Appendix). We obscured the spectacle scales of the snakes using nontoxic face paint (Snazaroo Classic; SI Appendix) to focus on control modalities where the animal reacted to collisions rather than avoided them. This behavior is likely relevant to this species as they do not appear to rely on vision during fast movement. This also prevented reaction of the animal to cues external to the experiment like movement of the researchers.
We focused our analysis on trials in which the snake initially traveled parallel to ẑ , passed through the array, then continued along a straight trajectory (Fig. 2 B and C). The pattern of alternating body bends was preserved throughout (Fig. 2 C and D). We measured an average decrease in speed of 84 ± 27% from before initial contact to once the snake was entirely clear of the array (SI Appendix, Fig. S4). While significant (P < 0.001, Wilcoxon signed rank test), the decrease was similar to that measured when no posts were present (83 ± 19% P < 0.001) and the animals were never observed stopping or turning back after contact.
After transiting the array, many snakes were deflected from their original heading (Fig. 3A). We characterized the pattern of reorientations by finding the angle with respect to ẑ , θ, of each trajectory (Fig. 2C) and calculating the probability density of all trials combined.
When no posts were present, θ was at most 25 • (Fig. 3A, Inset). In contrast, trajectories which passed through the posts were spread over -57.2 • to 56.1 • ; the animals were diffracted by the interaction like waves passing through narrow apertures (Fig. 3B). We observed three central peaks in the distribution, reminiscent of the interference pattern observed in matter-wave diffraction. These features were present in 80% of bootstrapping-predicted distributions (Materials and Methods) and were robust to the measurement method (SI Appendix, Fig. S5 B and D).
Reaction forces from the posts, Fpost (Fig. 4 A and B), persisted over timescales comparable to the time it took C. occipitalis to transit the array (Fig. 4 B and C; P = 0.29, signed rank test comparing 70% of the snake transit time to the total time forces above 3 mN were measured; SI Appendix, Fig. S7D).
We predicted the forces experienced by C. occipitalis body segments during movement across sand using resistive force theory (RFT), previously used to predict diverse locomotion in granular media (30) (SI Appendix). The median force measured by the posts was 5.2 times that calculated from RFT, which is similar to the ratio, 4.9, of drag force in the sandmimic versus sand (SI Appendix, Fig. S7 B and C). Thus, forces from the post were nonnegligible compared with those used for locomotion, implying that the snake did not change its selfdeformation pattern upon detecting a collision to avoid the obstacle.
Measurement of the force orientation angle, θ F (Fig. 4A), revealed Fpost was most often oriented perpendicular to ẑ , that is, the snake pushed "left/right" rather than "forward/back" (Fig. 4D). As the snakes were not observed changing kinematics to "grab" the posts or applying forces antialigned with the direction of motion at contact (i.e., θF = ±180), we concluded that the animals did not alter strategy upon collision to propel themselves forward using primarily the post. Interestingly, RFT predicted a similar distribution of θF arising from interaction with a sand substrate (Fig. 4D). While further study of the forces acting between the animal and the sand-mimic substrate would be necessary to make a definitive statement, this pattern of transverse forces is consistent with noninertial undulatory propulsion, suggesting that the animal interacts with the posts in the same way it pushes against piles of sand.
We hypothesized that forcing from the posts caused passive shape changes which broke the symmetry of the forward wave, causing reorientation of the snake (31). However, changes to the waveform during interaction with the post array were subtle (e.g., Movie S2 and Fig. 2 C and D). Therefore, to understand how the body deformed during interactions with obstacles, we elicited large deformation by replacing the post array with a solid wall. Vision was not obscured during these trials and snakes made no observable effort to avoid the collision, providing further evidence that these snakes do not rely on visual path planning during movement.
Upon collision with the wall the body "buckled" at the nearest extrema until the head traveled parallel to the board, at which point the animal changed strategy to bypass the obstacle (Movie S3 and Fig. 5A). This pulse of high curvature (Fig. 5B) was spatially localized; κ m of the next extrema, measured simultaneously, was unchanged (Fig. 5C). 2 of the central post. We calculated θ for each trial using data in the arc. (Inset) Probability density of θ measured when no posts were present. Trials were shifted such that the average direction of motion of the first third of all trials aligned with ẑ. The last two-thirds were used to calculate θ as in the array trials (n = 44, N = 8). The solid curve is the normalized probability density (integral under the curve is one) of θ. The shaded area is the 15.9th and 84.1th quantile estimated from a 10,000-iteration bootstrapping. We note that the small features in this curve were not significant (NS); they were less than the SD estimated using bootstrapping and were only present in 25% of the bootstrapping-predicted distributions (Materials and Methods), and their appearance depended on the measurement method and bin size used (SI Appendix, Fig. S5 A and B). BIOPHYSICS AND COMPUTATIONAL BIOLOGY PHYSIOLOGY C F ( °) B D (mN) time (s) F 1 F 2 0 1 20 0 -20 z x RFT experiment 0.015 -180 0 180 -90 90 0 0 0.5 1 time(|F|>3mN)/time(transit) 0 2.5 Probability Density A x z Fig. 4. Pattern of forces during transit through array. (A) θ F is the angle between Fpost and ẑ. Reaction forces F 1,2 were measured using post deflections (Materials and Methods). Diagram is not to scale. (B) Example post forces in x (top) and ẑ (bottom). The animal in this trial contacted a post on both its right (orange, F 1 ) and left (red, F 2 ). (C) Compare the force contact time-the total time in a trial the force magnitude |F| was above a threshold of 3 mN-to the transit time-the time between when the first tracked point first reaches the array to the last tracked point exiting (n = 233, N = 8). A value of 1 indicates the snake experienced force ≥ 3 mN for the duration of the time it was passing the posts. (D) θ F measured in experiment (black) and predicted by RFT (gray) (n = 194, N = 8, 28497 total measurements).
Undulating and Buckling Model. We previously observed "mechanical diffraction" patterns in an open-loop snake robot which was rigidly rotated to certain θ during interaction with posts (32). We hypothesized that mechanical diffraction in the animals could arise from adherence to a serpenoid self-deformation pattern which passively buckled in response to external forcing by the posts.
We developed a model to explore the outcome of our hypothesized strategy. Previous models of elongate, undulatory animals calculated the acceleration of body segments using internal forces, such as muscle activation and viscosity of the viscera, and external forces, like friction of the scales, discrete posts, or fluids (19,33,34). Rather than use a dynamical model of our system, we chose to develop a purely geometric model which allowed us to focus on the contribution and consequences of the hypothesized strategy.
We bypassed the complexity of the interaction between the spatially extended body, the substrate, and the post by assuming forces between the body and the substrate were such that the snake propagated with no slip and the body was buckled by the post the minimum amount necessary to prevent overlap. The position of each segment was prescribed at each time based on these assumptions (35)
During terrestrial lateral undulation snakes self-deform using waves of unilateral activation of the epaxial muscles (36). Thus, we assumed passive body buckling would occur at those locations on the waveform where the shape change would further shorten active muscle segments (37). We dictated the model κ changed at a "preferred buckling" location, s buckle , determined by this muscle activation pattern (Fig. 6A).
When the model snake contacted a post, curvature at s buckle was set to the predetermined absolute maximum value of κ abs = 50 m -1 (nominal amplitude 25 m -1 ) over as many segments necessary to solve the constraint (Movie S4), leading to a pulse of high curvature like that observed in the snake wall-collision trials (Fig. 6B, Inset). Given the assumption that the body was passively deformed by external forces, the modified waveform was passed down the body at the wave speed. In both the model the biological snake, interactions with each post were independent.
We initiated the model trajectories at 700 locations encompassing all possible interactions (SI Appendix, Fig. S8A). Like C. occipitalis, the model trajectories were spread by the array interaction (Fig. 6B), either continuing approximately along ẑ or at an angle θ = ±22.0 ± 7.0 • (Fig. 6C). There was a relationship between the initial conditions and the scattering angle, suggesting this system is deterministic (SI Appendix, Fig. S8D).
The location and relative prominence of the peaks depended on the geometric parameters. We defined R buckle = (κ abs rpost) -1 and D buckle = smaxr -1 post , where smax = L 2ξ was the largest possible value of s buckle . In the model we used the same d = 23 mm as in the snake experiment but increased post radius to rpost = 10 mm as the model was infinitesimally thin. We used R buckle = 4.1 and D buckle = 20 in the model. These parameters were comparable to the average snake values of R buckle = 3.5 ± 0.4 [mean and range estimated using both the wall trials and an anesthetized snake (38)] and D buckle = 31.8 ± 3.1. This combination of parameters yielded the same distribution of post contacts as in experiment (Fig. 6C, Inset and SI Appendix, Fig. S10A).
Inspired by the effect of slit width on diffraction patterns, we used the purely deterministic model to examine θq , the spread of the distribution, as R buckle and d varied. We found as R buckle increased θq decreased (Fig. 6D). Increasing R buckle meant that more segments were involved in buckling to bypass the post, effectively increasing s buckle and requiring less angular deflection of the trajectory (as seen in Fig. 7B).
Similar to diffraction of fluids or subatomic particles, θq decreased as d increased (Fig. 6E). When freely moving, the maximum angle the body made relative to ẑ was φo = 45.6 • . At d = 23 mm the maximum angle the body could make between the posts was φp = 29.6 • and θ was determined primarily by the wave phase at contact (SI Appendix, Fig. S8D, Top). In contrast, when d = 50 mm φp = 66 • , the model never contacted more than one post in a trial, and the impact location on the post was influential (SI Appendix, Fig. S8D, Bottom). We rationalized that for small spacings there was a greater probability the body had to reorient to match φp regardless of the initial contact location. As φp became larger than φo the wave always "fit between" the resists lengthening does not resist lengthening does not resist shortening post post ON ON OFF 1v o T 0 100 pdf v o T posts A B C D posts and θ were dominated by the trajectory deflection around the initially contacted post. Further, as d increased more trials transited the posts without contact, increasing the signal at θ = 0 and moving θq inward.
Ray Model. To gain insight into the fundamental interactions governing snake reorientation we developed a simplified analytical description: we represented the traveling wave as a "ray"; the limiting case of the serpenoid curve (sinusoidal κ) as κm → 0. The ray interacted with the post following rules on buckling locations and the bending limit derived from the buckling model and experimental measurements, respectively. Such a reduced approach is beneficial in simplifying calculations for scattering in wave systems [e.g., specular scattering of light (39)].
For the current study, we focused on understanding interactions with a single post. While we were not able to perform the snake experiments for a single post due to the low probability of the animal's contacting an obstacle, using the undulating model we noted that the width and peak locations in the distribution of θ were comparable to the multipost case (SI Appendix, Fig. S9C). The most noticeable difference was the higher central peak relative to the side lobes, likely because the addition of multiple posts remapped trajectories that would have had small-scatter collisions with the single post to collisions which resulted in larger angles (32).
We calculated the distribution of scattering angles for rays initially parallel to ẑ which bent to κ abs at location s buckle to solve the post constraint (Fig. 7A). The predicted distribution agreed well with that calculated using the undulating model for a single post (Fig. 7B). Scattering angles from head-buckle (s buckle = 0, Fig. 7A, Left) trajectories rose from zero probability at θ = 0 to a maximum value set by θmax = cos -1 ( κ -1 abs (κ -1 abs +rpost)
) before dropping off, yielding the secondary peaks. The central peak was primarily trajectories which buckled at locations s buckle > 0 (Fig. 7A, Right).
The agreement between the ray model and the more complicated undulating model suggested that the features of the scattering pattern were more dependent on the physiological constraints (bending limit and a spatially extended body) than the periodic waveform. Future work could include extending the ray model to include multiple posts.
Our surprising discovery of a mechanical diffraction pattern arising from interaction between a sand-specialist snake and multicomponent terrain revealed the benefit of passive mechanics in negotiating collisions in limbless systems. We used a geometric model to show that the reorientation pattern was reproduced by an open-loop control strategy wherein the motor program continued unaltered during unexpected collision with the posts while unilateral muscle activation allowed the body to passively buckle around obstacles. A ray model provided insight into the fundamental mechanisms behind the diffraction, implicating that the existence of preferred buckling locations on an extended body, engendered in our case by the unilateral muscle activity pattern, caused trajectories to scatter nonuniformly off of a post.
Our work provides a starting point for disentangling the role of mechanics and active control in limbless terrestrial systems; testing hypotheses for neuromechanical control in a frame as suggested by ref. 40 would be illuminating in this regard. The use of mechanical structure to supplement open-loop control has been largely studied in rapidly moving legged animals (9,14) and used to simplify control in bipedal and sprawled-posture robots (41,42). This work suggests that passive mechanisms are quite general, and such strategies could be useful in the next generation of limbless robots where the control paradigm has
been closed-loop tracking of joint trajectories (43). However, future work is needed to determine when this strategy is not appropriate (e.g., terrains with more dense obstacles or even certain waveform/obstacle combinations). Finally, we comment on our results in the context of the emerging field of active matter and collisions which do not conserve momentum (44). Our biological mechanical diffraction mimics phenomena found in subatomic systems; this is part of a growing realization that active collisions are a fertile source of interesting dynamics from mechanical diffraction in a robot (32) and scattering of microorganisms off walls (45,46) to shape-induced reorientation of a cockroach (8). We posit that a framework which takes inspiration and tools from diverse systems will broaden understanding of principles of heterogeneous interaction in self-propelled systems across scales.
All codes and datasets are available at smartech.gatech.edu. Snake Experiments. All C. occipitalis were collected in accordance with scientific collection permits (nos. SP790952, SP625775, and SP666119) approved by the Arizona Game and Fish Department. All snake experiments were conducted under the Georgia Institute of Technology IACUC protocols A14066 and A14067. Snakes were set in the arena and once they began moving were not contacted during a trial. Temperature in the testing and holding area was maintained at 27.2 ± 0.6 • C. All reported values are mean± SD unless otherwise indicated.
Bootstrapping. We used bootstrapping to estimate the importance of the features observed in the θ histograms. Using θ measured in 194 trials we generated 10,000 distributions using random resampling with replacement. We counted how many times the three central peaks and two central valleys were present in the distribution by searching for features of prominence at least 0.004 (smallest prominence in the real distribution = 0.0052) occurring on the interval between its neighboring features. Peak locations, determined by their locations in the distribution were [-30, 10) • , (-10, 15) • , and (15,35] • and valleys were (-20, 5) • and (5, 20) • . The same procedure was used to count the three peaks in the distribution without posts with peak locations at [-30, -5) • , (-5, 10) • , and (10,35] • and valleys were (-10, 5) • and (5, 15) • .
Force-Sensitive Posts. The theoretical prediction for force, F, as a function of peg tip deflection, δ, is F =
3Eπr 4 post 2a 3 (9 L post a -5)
δ. Lpost = 7 cm, a = 2.81 cm is the height of applied force measured from base, Young's modulus E = 5.7 ± 0.6 MPa. Peg tip deflection δ was measured from the high-speed video. Post deformations were small compared with snake length scales (3-mm deflection at tip for a load of 0.050 N at a = 2.7 cm). Uncertainty due to variation in a was < 10% (SI Appendix, Fig. S3C). Relative movement between the arena and the camera was corrected using the fiducial post.