Mitigating memory effects during undulatory locomotion on hysteretic materials
Abstract
While terrestrial locomotors often contend with permanently deformable substrates like sand, soil, and mud, principles of motion on such materials are lacking. We study the desert-specialist shovel-nosed snake traversing a model sand and find body inertia is negligible despite rapid transit and speed dependent granular reaction forces. New surface resistive force theory (RFT) calculation reveals how wave shape in these snakes minimizes material memory effects and optimizes escape performance given physiological power limitations. RFT explains the morphology and waveform-dependent performance of a diversity of non-sand-specialist snakes but overestimates the capability of those snakes which suffer high lateral slipping of the body. Robophysical experiments recapitulate aspects of these failure-prone snakes and elucidate how re-encountering previously deformed material hinders performance. This study reveals how memory effects stymied the locomotion of a diversity of snakes in our previous studies (Marvi et al., 2014) and indicates avenues to improve all-terrain robots.
Introduction
Movement is crucial to the survival of most organisms and a necessary ability in robots used in fields like medicine (Taylor, 2006), search and rescue (Murphy et al., 2008), and extraterrestrial exploration (Lindemann and Voorhees, 2005; Shrivastava et al., 2020). Successful locomotion depends on the execution of body-shape-changes which generate appropriate reaction forces from the terrain, a relationship that can be further complicated if the body motion permanently changes the state of the substrate. Much of our knowledge of terrestrial locomotion is in the regime of rigid materials where the terrain is not affected by passage of the animal (Koditschek and Full, 1999; Sponberg and Full, 2008; Kelly et al., 1997; Farley et al., 1993), an understanding which led to the development of robots which are effective on hard ground (Saranli et al., 2001; Liljebäck et al., 2012; Blickhan et al., 2007).
Little is known about locomotion in non-rigid materials which are plastically deformed by the movement of animals or robots, leaving tracks or footprints. Deformable terrains represent a spectrum from materials like water which continuously flow toward the undisturbed, zero shear state (Lautrup, 2011) to those which are remodeled by the interaction, like sand. Previous work has mainly focused on motion in fluids, where disturbances dissipate (Ijspeert, 2008; Sfakiotakis et al., 1999; Liao, 2007); the impact of soft material hysteresis on locomotion is not well-understood (Mazouchova et al., 2013; Zhang and Goldman, 2014).
At one end of the spectrum, where deformations are short-lived (Boyle, 2010), small fluid swimmers like the nematode Caenorhabditis elegans (Wen et al., 2012), spermatazoa (Gray, 1953), and bacteria in water (Rodenborn et al., 2013) are propelled by the viscous force of the material resisting the animal’s body shape changes. In these low-Reynolds-number systems (Re, the ratio of inertial to viscous forces), the inertia of the microscopic animal is negligible compared to the fluid viscosity such that if the animal stops self-deforming it very rapidly stops translating and/or rotating. Such resistive-force-dominated swimming is also found at the macroscale in frictional-fluid swimmers like the sandfish lizard Scincus scinus and the shovel-nosed snake Chionactis occipitalis moving subsurface through granular matter (GM) (Maladen et al., 2009; Sharpe et al., 2015). The motion of these animals is similarly dominated by resistive forces, in this case the frictional interactions between grains. Notably, in both viscous and frictional fluids, the material surrounding the submerged swimmers continuously re-flows around the body of the animal such that the animal is always surrounded by material.
At the other end of the spectrum, animals move in materials whose state depends on the history of interaction. Swimmers in fluids at high-Re, where inertia dominates fluid viscosity, can experience time-dependent flow (Smits, 2019; Oza et al., 2019). Less well-understood is how animals manage material memory (Keim et al., 2019) when traversing terrestrial substrates like soil, mud, or sand. Such materials will yield in response to forces applied by a locomotor, similar to a fluid. However, unlike fluids, soft materials can typically bear internal stress. At the free surface, gravity is often insufficient to return mounds of material created by a locomotor to the undisturbed state, exemplified by the tracks left behind by a passing animal. Our previous work (Mazouchova et al., 2013; McInroe et al., 2016) revealed that performance of limbed animals and robots is affected by the ability to avoid interactions with previously-disturbed material. Using limbless lateral undulation on deformable substrates to generate appropriate propulsive forces is a non-trivial task, as demonstrated by the failure of both a snake-like robot (Maladen et al., 2011a) and a number of snake species (Marvi et al., 2014) attempting to traverse the surface of GM. However, it is unclear how such animals (or robot models) might manage environmental interactions to move effectively when, unlike limbed locomotors, they cannot increase step size to avoid their own tracks.
In this paper we gain understanding of movement through complex terrestrial materials by combining several approaches, each of which provided insight into the number of interrelated elements governing the system. We begin with a study of a variety of non-sand-specialist snake species moving on the surface of GM as well as a desert-specialist snake. We quantify the kinematics of the snakes, which move with varying degrees of success using a range of waveforms. To understand the connection between waveform and performance, as mediated by the GM, we measure the granular response to surface drag. We develop a model for speed dependence that indicates that the granular physics is unchanged over the range of speeds observed in the animal. We then find a characteristic drag anisotropy curve that is largely independent of drag depth and speed. Based on the granular drag experiments and our animal observations we hypothesize that, despite the fast movement and complex granular flows, body inertia is negligible. Thus, we introduce a surface granular resistive force theory (RFT). Using RFT, we identify a trade off between actuation speed and torque from the GM that explains the stereotyped waveform used by the sand-specialist as that which maximizes movement speed under an anatomical power constraint. We show that performance of the non-sand-specialist snakes depends on both morphology and waveform, although RFT only accurately predicts performance when lateral slipping of the body is less than a threshold value. In these cases we often observe the animals re-interacting with material remodeled in previous undulation cycles. We study the interplay of waveform and material remodeling using a robophysical model. RFT calculation of robot locomotion, as in the high-slip snakes, is inaccurate. Particle image velocimetry (PIV) measurements indicate RFT over-estimates performance when the robot body is pushing material lateral to the direction of motion. RFT does, however, provide insight into the regime of locomotor failure, where we find the desert-specialist operates far from failure due to a combination of its waveform, slender body, low-friction scales, and lifting of segments that only produce drag. We close with a summary of our results and a brief discussion of the implications.
Better understanding of the connection between body self-deformations, substrate remodeling and force generation, and the resulting locomotor performance will elucidate principles for effective motion on hysteretic materials. This can be applied to the next generation of all-terrain robots. Further, future studies can leverage such knowledge of the benefits and limitations of different body–terrain interaction modes to tease apart neuromechanical control strategies for contending with natural terrains (Schiebel et al., 2019).
Results and discussion
Performance of limbless locomotors on granular material
Snakes occupy a variety of habitats and display a wide range of morphologies. We took advantage of this natural diversity to explore how different body shapes and patterns of self-deformation fare on GM (granular matter, see Table 1 for a list of symbols and abbreviations). We studied 22 species of snakes representing five families from the collection at Zoo Atlanta. Snake body plans ranged from short and stout to long and slender (Figure 1(a,b); Table 2), and their natural habitats encompassed a broad range from wetlands and swamps, to wet and dry forests and rainforest canopies, to rocky deserts and mountains.
As a counterpoint to the variety of snakes, which were either terrain generalists or specialized to habitats which did not have an omnipresent granular substrate, we also studied the shovel-nosed snake Chionactis occipitalis (Figure 1(c); Table 3). This species is specialized to bury within (Kavanau and Kavanau, 1966; Sharpe et al., 2015) and move across (Mosauer, 1933) the dry sand of their desert habitat. We used nine individuals collected from the desert near Yuma, Arizona, USA (see Materials and methods for details). C. occipitalis uses a stereotyped waveform (Schiebel et al., 2019) to quickly move across the sand with little slipping of the body (Mosauer, 1933).
Kinematics and ability vary among species
The non-desert specialist snakes from the collection at Zoo Atlanta were tested in a 2 × 1 m2 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 cm2 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, , 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),
The speed that the wave propagates down the body, , determines how quickly the shape changes in time, . At any moment, the shape of the snake as a function of arclength along the body, , is determined by the amplitude, , also referred to as the attack angle, and spatial frequency, or number of complete waves on the body, (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 and 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 and . 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 and representative of the average snake posture. Plotting the results for each of the nine C. occipitalis in the versus 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 and when considering both each individual and each trial (Figure 2(b), inset). Given the weak relationship between and , 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.
Granular drag measurements
Principles governing subsurface swimming in GM were elucidated by granular drag measurements (Maladen et al., 2009). These experiments revealed that during subsurface swimming in GM the material behaved as a frictional fluid. Like a low-Re swimmer in a viscous fluid, propulsive forces arose from the resistance of the GM to the animal’s shape changes (although rather than viscous forces the grains interact with each other via normal and frictional contacts Maladen et al., 2009; Zhang and Goldman, 2014).
We observed the formation of granular piles as the snakes self-deformed on the granular surface (Figure 3 (a), Figure 3—video 1). This suggested that, consistent with subsurface swimming (Maladen et al., 2009) and surface walking (Mazouchova et al., 2013; McInroe et al., 2016), slithering animals propel themselves using the granular forces arising from the body yielding the material (and not solely frictional anisotropy of the ventral scutes Hu et al., 2009). The interaction between the body and the GM can be characterized by the amount of slip, , relative to the substrate, where and are the local velocity and tangent unit vectors of a body segment (Figure 3 (a) inset; Sharpe et al., 2015). emerges from the interrelationship between the granular stress and the animal’s self-deformation pattern.
We empirically measured the granular stress on a simple model for a snake body segment—an aluminum plate commanded to move at drag angle between the drag velocity unit vector and the plate face tangent in the horizontal plane (Figure 3 (a,b)). An aluminum rod attached the Al plate to a force transducer (ATI nano43) that decomposed stresses into those acting normal , and tangent, , to the plate face (Figure 3 (c)). This apparatus was located on the end effector of a six degree-of-freedom robot arm (Denso VS087A2-AV6-NNN-NNN) which controlled intruder movement. The robot arm first rotated the plate to , then submerged it to depth , measured from the free surface to the plate’s bottom edge (Figure 3 (b)), next it dragged the plate parallel to the surface for 20 cm at constant speed, and lastly stopped and extracted it. A fluidizing bed containing the same 297 ± 40 μm glass particles used in the C. occipitalis experiments prepared the material to an initially loose-packed state using two shop vacs controlled by a proportional relay. The fluidization procedure was the same as in the snake trials.
Unlike subsurface drag stresses, which developed almost instantaneously to the steady state (Maladen et al., 2009), at the surface stress monotonically increased over several centimeters before saturating (Figure 3(d)). This is due to the free surface flow of the GM; a pile of sand above the surface, like those created by the snakes, appeared at the leading face of the intruder at the onset of drag and increased in volume until reaching a balance between the new grains being encountered and those flowing around the edges of the plate. Previous studies of plate drag at the surface (at deg) measured a similar drag-distance dependent force and observed a wedge-shaped region of GM, beginning at the bottom edge of the plate and extending above the surface, whose constituent grains were flowing forward and up against gravity (Gravish et al., 2014; Guo et al., 2012).
Velocity-dependent stress captured by grain inertia model
The snake speeds were variable ( 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 deg and mm and found normal stress quadratically increased as increased from 1 mm s-1 to 750 mm s-1 (the limit of the robot arm capability, Figure 4(a)). Similarly, for deg and mm s-1, normal stress linearly increased with 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)).
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 term in the formulation is the momentum transfer requirement for pushing the granular material ahead of the plate element, and the static () stress term, , is the stress required to initiate yielding of the media. is the density of the material, which we estimated using the density of glass, and the packing fraction induced by motion of the intruder, .
The static stress on the intruder may be modeled using gravitational loading of the grains. In fluids, hydrostatic pressure is given by given gravity and local plate element depth . In granular materials there is a similar depth dependence, however, because the bulk can harbor internal stresses, the behavior differs substantially and where 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 , so . We assume is to estimate that at mm ( mm) Ncm-2 and at mm Ncm-2 which are similar to the experimentally measured values.
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 estimated from the experimental data. was calculated by subtracting from measured at the three lowest speeds collected at a given and averaging the result.
Using the empirically obtained static stress, the term captured stress as a function of 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, , 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 (Figure 4—figure supplement 1(a)) while the waveform was independent of (Figure 4—figure supplement 1 (b,c)). This suggested that, despite the 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.
Drag anisotropy is not strongly dependent on speed or depth
At low-Re, the normal and tangential stresses acting on segments of long, slender swimmers in Newtonian fluids are approximated by and , respectively (Gray, 1955). The constants and are determined by the geometry of a body segment and the viscosity of the surrounding fluid, and the ratio 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 (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, . 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 ranging from 0 deg to 90 deg (Figure 3(b)). The ratio was largely independent of the drag distance (Figure 5(a)), especially at small . There was a periodic fluctuation of occurring over several cm appearing in all trials. However, both these fluctuations and the slope of as a function of drag distance were small relative to the effect of changing . Thus, we averaged and to characterize the relationship between plate orientation and stress normal and tangential to its face (Figure 5(b)).
Normal stress increased monotonically with while tangential stress was approximately constant, gradually falling to zero as went to 90 deg (Figure 5(b)). The ratio as a function of 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.
In viscous fluids the constants and are independent of drag angle such that is constant. Consistent with results for subsurface drag (Maladen et al., 2009; Maladen et al., 2011b), we find that at the surface is a nonlinear function of (Figure 6(b)). Because tangential stress magnitudes were relatively small and, unlike in viscous fluid, normal stresses rose sharply with at small angles (Figure 5(b)), thrust and drag forces on the plate were equal () at smaller than in viscous fluid or subsurface in GM (Figure 6(a)). Reflecting the small angles where the anisotropy was unity, measured on the snake were small (Figure 6(a), gray histogram).
Plotting anisotropy at constant as speed and depth varied further illustrated the relatively weak dependence of on these variables (Figure 6(c,d)). Especially at small , anisotropy did not depend on (Figure 6(c). was more dependent on at shallow depths, decreasing twofold as increased from 4 to 8 mm (Figure 6(d)).
The invariance of 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 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.
Surface resistive force theory model
Continuum constitutive equations for granular media (similar to Navier-Stokes equations for fluids) can predict granular flows (Kamrin and Koval, 2012). However, these methods, as well as established discrete element methods (like MD), are computationally expensive. Granular resistive force theory, which relies on the assumption that the total force on the body is a linear, independent superposition of the forces on individual segments (Gray, 1955), is by comparison computationally cheap. RFT successfully modeled a number of systems for which these assumptions are valid, although its effectiveness in a system exhibiting hysteresis was unknown (Zhang and Goldman, 2014).
We observed that as a trunk segment pushed laterally against the sand it built a pile of grains, like that which evolved during drag of the plate, on the side of the body closest to the tail (Figure 3—video 1). The side of the body opposite the pile (closest to the head) appeared to be in contact with little material. Therefore, we chose to model the snake body as a flat plate, representing the caudal-facing side of the body, with an added term to account for drag on the ventral surface. This is not an unreasonable model for this species, as the ventral surface is known to be flat or even slightly concave (Mosauer, 1933). As the snake held the head slightly raised from the surface we did not include a term accounting for head drag.
We approximated the kinematics of the snake using segments evenly distributed along the arclength () at 70 points in time () divided evenly over one full undulation cycle. At this resolution adding or subtracting five segments or time steps changed the velocity prediction by less than 0.1%. We calculated each segment’s orientation and velocity given parameters and using the equation for a serpenoid curve (Equation 1). The segment speed, , was 1 ms-1 unless otherwise stated.
The force on a segment of length was assumed to depend only on its movement through the material, , given by its orientation (, ) and direction of motion, (). Velocity, , was determined by the local segment velocity, , added to the center-of-mass velocity, , and rotation about the CoM, . and were calculated using Equation 1. Given functions and relating granular force normal and tangential to the segment, respectively, to the segment’s motion, the infinitesimal force was thus
z was the plate intrusion depth (as defined in Figure 3), was the width of the plate used in the granular drag experiments, and was the snake intrusion depth. Unless otherwise noted, we assumed mm. This prefactor scales the granular forces under the assumption that force was linearly related to intruder area. While shape effects on drag at the surface are not well known, during subsurface drag the horizontal forces on an object depend primarily on cross-sectional area with little dependence on the shape of the object, much less than in fluids (Albert et al., 2001). There are also shape dependent lift forces on an immersed intruder (Ding et al., 2011), however, we assumed that vertical forces acting on the snake were balanced and did not affect the horizontal forces determining performance. Our depth-dependent drag experiments, which revealed drag stress increased linearly with depth (Figure 4), also indicated force scaled linearly with area.
The function was given by
where , and were constants determined by fitting the drag data. For example, at mm and mms-1, (i.e. data in Figure 5(b)) N, , and . This is similar to the function presented in Sharpe et al., 2015, with the addition of the overall scaling constant . We found that was not well-captured by our previous subsurface models. In the interest of carrying out the RFT calculation, we thus used a Fourier fit to the data to determine (Table 4).
We modeled ventral drag on each segment as acting opposite the velocity with magnitude where μ was the coefficient of friction between the snake’s ventral scutes (scales on the belly) and the substrate, was the snake mass, and was the gravitational constant.
Previous research (Hu et al., 2009) found that the ventral scutes were anisotropic such that gliding directly forward produced less drag than sliding laterally, and both forward and lateral motion resulted in less frictional drag than moving directly backwards. The coefficient of static friction between C. occipitalis scales and the glass beads was found to be 0.109 ± 0.016 when sliding forward and 0.137 ± 0.018 backward (Sharpe et al., 2015). The lateral coefficient of friction for C. occipitalis scales is not known. However, given the results of Hu et al., 2009, we assumed that it was bounded between 0.10 and 0.14. Thus, the difference between frictional forces acting tangential and normal to segments was negligible compared to the forces due to the motion of the body through the GM Therefore, we chose not to include the frictional anisotropy in our RFT calculation.
We previously found the coefficient of static friction between aluminum and the glass particles was approximately 0.2 (Maladen et al., 2009). To account for the difference between the plate friction and snake-scale friction, we scaled the tangential drag forces by the ratio of the scale friction to the plate-GM coefficient. For example, when predicting the performance of C. occipitalis, we scaled the function measured using the aluminum plate by 1/2.
The RFT calculation was carried out as follows. At each time step and of each of the 100 were determined using Equation 1. Given some and , can be calculated from Equation 3 for each . Summing over the body thus gives the total force acting on the CoM in and and torque about the CoM. Assuming the animal was in the steady state, such that the total force/torque acting on the body was zero, MATLAB’s lsqnonlin function was used to find the and components of and (assuming planar motion) component of which resulted in zero net force and torque at each of the 70 time steps.
Surface granular RFT accurately predicts sand-specialist speed
We began by using RFT to estimate the steady-state of the sand-specialist given the average ms-1 measured in experiment, snake length 39 cm, body radius 3.5 mm, and mass of 19 g. We calculated the average by integrating 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 (Figure 7(a)). Notably, we did not account for material hysteresis in our calculation.
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 to larger (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 near 0 deg, where (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 (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 and less than the experimental uncertainty in (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 , and captured the relationship between and 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.
Trade-off between internal and external factors constrains waveforms
While the total torque about the CoM, summed over all , at each time was zero, the torque on any given segment was not. We thus used RFT to predict the maximum joint torque, , experienced by any C. occipitalis body segment over an undulation cycle (Figure 8(a)).
We calculated torque acting on the CoM assuming infinitesimal width and finite segment length. For joint the internal torque was (Ding et al., 2011). We find this value for each of the segments at all times. is then the maximum value in this 100 × 70 matrix.
increased as or decreased. This was both because balancing forces required larger , 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 and/or , 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 : increase the wavelength by decreasing and/or (see example waveshapes in Figure 8(b)) or increase the speed of wave propagation (given the linear relationship between propagation speed and , 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 , we explored how quickly a ‘muscle’ segment of fixed nominal length relative to total body length would have to shorten as a function of and . 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 for that shape’s stride length. We modeled a body segment as a bending beam with length along the spine and arclength along the inside of the curve (see diagram Figure 8(b)). We approximated the amount of shortening the muscles must produce as the difference, , at the point of maximum bending. Using this model we estimated , the speed the inside of the bend must change length to go from at the inflection points of the wave to the maximum occurring at the peaks given undulation period 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, (which we take to be that occurring at the mean and measured in experiment). The actuator speed increased with and . 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 and/or 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.
C. occipitalis waveform maximizes movement speed under anatomical constraint
We endeavored to include the trade off between decreasing torque and increasing actuation speed needed to maintain (now using RFT to estimate the actual stride length given slipping) as and increased. To do so, we used RFT to calculate the joint-level power of each segment over a cycle for each shape in the , space.
The angle between adjacent segments, , relates to via . Joint-level power for each joint at each time was calculated using the rate of change of the joint angle, , where
The power-limited velocity, , 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, was maximized in the center of the , space inhabited by the waveform of C. occipitalis (Figure 9(a)). Reflecting the oscillations in , this metric had maxima near integer values of (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 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 reflects differences in peak muscle power capabilities.
Stout snakes must use larger attack angles to succeed
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, , 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, (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, , and 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, kgm-3, and a circular cross section. The width and length of a snake can then be calculated as and .
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 (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).
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 the snakes were, with the exclusion of one exception, moving effectively with an average deg (Figure 10(a), blue crosses). Below 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 (Figure 10(a), red crosses). The other snakes of moved with higher slip on average than their high-aspect ratio counterparts, deg (Figure 10(a), gray crosses).
We found that animals with 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 changed for those snakes with an aspect ratio less than 26 (and the snake that failed with =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 N = 6, n = 4311 measurements, N = 12, n = 9017, p << 0.01). Those snakes which had 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 deg and 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 and slip 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 ( 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.
Material remodeling regulates locomotion
The waveform and performance of the snakes was variable, and this variability was reflected in the tracks left by the animals. As the material did not re-flow around the body after being disturbed, we hypothesized that the manner in which different waveforms remodeled the substrate was important in determining performance. We systematically explored the impact of waveform on performance using a robophysical model, a 10 joint robot on the surface of poppy seeds (Figure 11).
Robophysical model elucidates connection between waveform and performance
The robot executed serpenoid curves (Equation 1), realized using a Lynxmotion SSC-32 servo controller to dictate the angles between adjacent motors, , (Figure 11(a)). Each motor was offset from its anterior neighbor by a constant phase set by . We varied and and recorded kinematics as the robot performed three undulation cycles on the surface of poppy seeds in a 1.2 × 1.8 m2 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 between 1 and 1.4 because lower 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; 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.
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 and (e.g. Figure 9), and was commensurate with our findings in the biological experiments (Figure 10 (b)). Over the range of available to the robot, performance was more strongly dependent on attack angle than wavenumber, therefore we focused our attention on variations in .
RFT prediction of robot performance is inaccurate
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 for fixed (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 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.
Material remodeling changes 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 , 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 over a run for and each of the three tested (Figure 13(b)).
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 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 waveforms had a significant velocity component lateral to the direction of motion; as decreased, the peaks in the 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 , 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.
RFT predicts performance transition
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 , space which RFT predicted to be at deg separated the successful and unsuccessful robot waveforms (Figure 14(a)).
We found there were a limited range of waveforms predicted to both have 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 and ). 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.
Conclusions
Undulatory motion on the surface of materials with memory is challenging for both animals and robots because the sensitive connection between a locomotor's waveform and its performance is complicated by material remodeling. We begin the search for principles of locomotion on substrates with memory by studying a sand-specialist snake, generalist snakes, and a robophysical model. We also used new granular drag measurements to extend resistive force theory to calculate surface performance. By comparing the limbless locomotors and using surface RFT calculations, we identified a number of factors which improved undulatory locomotion on our hysteretic material—sufficiently high attack angle waveform, high aspect ratio body, low scale friction, and lifting of the wave apexes. The generalist snake and robot experiments demonstrated that possessing all of these is not necessary for motion. However, as morphology becomes less optimal (more stout, higher friction) effective management of the substrate remodeling via choice of the appropriate waveform was increasingly crucial to progress.
The interrelationship between morphology, waveform, terrain remodeling, and performance rationalizes why it has been difficult to develop slithering robots which move on terrestrial materials with memory as compared to swimming in fluids or using rigid obstacles. Currently, robots are limited by the power density of the motors to low aspect ratio morphologies with a limited range of available waveforms such that they are in a regime where small changes in waveform can impact progress through the environment. The continued development of smaller, high-power-density motors will allow snake-like robots with higher aspect ratios and more joints. This study indicates that such robots will correspondingly be more adept at movement in substrates with memory, like those found in many natural terrains.
An immersed locomotor is aided by fluid-like re-flowing of the GM which mitigates memory-dependent failure modes (Maladen et al., 2011a). However, when C. occipitalis is on the surface its escape response is to flee across the GM rather than dive into it (Warren, 1953). We hypothesize this is because, despite the challenges we identified, there are a number of advantages for a locomotor on the surface. The granular stress magnitudes are less than subsurface and granular stress anisotropy rises more sharply such that the faster-moving, low attack angle shapes can be effective. Drag forces are reduced by avoiding the need to push the head through the material and the opportunity to lift segments which are not generating thrust. There is also less friction on the body, which at the surface is due only to animal mass instead of both mass and granular pressure. Worthy of note, locomotor performance is repeatable, suggesting the history-dependent forces are a predictable function of body shape changes and granular remodeling. Future work elucidating such a predictive model could both facilitate design of more effective all-terrain robots and further understanding of the biology.
Materials and methods
Animal experiments
Request a detailed protocolAll experimental procedures carried out at Zoo Atlanta were conducted in accordance with the Georgia Institute of Technology IACUC protocol A14058. All procedures involving zoo animals were reviewed and approved by the Zoo Atlanta Scientific Research Committee.
C. occipitalis were collected by Kevin and April Young in the Colorado Desert near Yuma, Arizona, USA under scientific collection permits (SP790952, SP625775, SP666119) approved by the Arizona Game and Fish Department and held in the Physiological Research Laboratory at Georgia Tech. Neither the sex nor the age of the animals was determined; gender and age dependent effects were beyond the scope of this study. All experimental procedures were conducted in accordance with the Georgia Institute of Technology IACUC protocols A14066 and A14067.
The temperature in the track way and snake holding area was measured prior to each trial. Lamps were used to ensure the temperature in both remained at 26 ± 1 °C , within the active range for C. occipitalis (Klauber, 1939). The heat lamps on the track way were turned off during data collection and LED lights were used for illumination.
Each day the C. occipitalis to be tested were transported from the housing facility to the lab where we conducted the trials. Snakes which were in the process of shedding were not used. During a trial, the snake was removed from its holding container and placed immediately in the fluidized track way. The snakes tended to be skittish and handling both during trials and in the housing facility was kept to a minimum. The animals would often immediately flee across the surface upon introduction to the track way; otherwise a light tail tap would elicit an escape response. If an individual did not respond to this stimulus they were returned to the holding container. If an individual failed to perform for three trials in a row they would be retired from the day’s studies. Snakes were tested at most every other day with a maximum of two successful trials collected per day. A run was included if the snake performed at least four complete undulations moving along a straight trajectory at apparently constant speed.
Calculation of waveform parameters
Request a detailed protocolTo measure we found the points of zero curvature; these corresponded to the inflection points of the waveform. We then calculated the arclength between these two points for the first half-wave on the body, multiplied by 2 to get the arclength of one full wave, , and divided the individual’s length by the result.
We calculated of C. occipitalis using finite differences to estimate the x and y values for each segment’s tangent vector. We generated sample serpenoid waveforms of known parameters and used these to determine that, in the presence of white Gaussian noise, the most accurate measurement of the tangent vector was obtained by subtracting the average position of the segment in question and three segments anterior from the average of the segment and the next three posterior segments. This method helped buffer against noise while still providing an accurate measurement. We found the maximum angle on the body at each time step and then averaged over all times in a trial in the to obtain .
Calculating of the non-sand-specialists was challenging as in some cases direction of travel was difficult to determine or the waveform was tortuous. Thus, we estimated using and the wave amplitude in curvature, , to determine the non-dimensional (Sharpe et al., 2015). Assuming serpenoid curves, . We used finite differences to calculate and the same procedure as described for to determine the best result was obtained by averaging over five points anterior and posterior to the point of inquiry.
Bending beam model for snake segment motion
Request a detailed protocolWe model a snake body segment as a beam with midline length and length of the inside of the curve . For a given body radius, , and radius of curvature, , . The speed the inside of the segment must change length from the unbent length of to is . We can write as . Using the relation we get . Next, using for a serpenoid curve , we find the maximum curvature in terms of attack angle . Thus the maximum length change . The time the segment has to undergo this change is a quarter of a period, so . Thus the shortening speed of the inside of a segment is
Lastly, for the nominal values we get the ratio and the relative shortening speed is
For a predetermined , MATLAB’s fsolve is used to find the undulation frequency for a given and which results in no-slip motion at that value of .
Dissection and torque estimation
Request a detailed protocolFour adult C. occipitalis (SVL 33.0 ± 3.7 cm, mass 16.4 ± 3.4 g) which had preserved in formalin and stored in ethanol were scalened and the body (from the posterior margin of the quadrates to the cloaca) was cut into 10 equal length segments, which were then weighed intact. Segments were then eviscerated and the epaxial muscle mass (consisting of the largest muscles, the m. multifidus, m. semispinalis-spinalis, m.longissimus dorsi, and m. iliocostalis Jayne, 1988) was removed, and the viscera, muscle, and remaining body tissue were weighed. All tissues were kept moist in 70% ethanol and dabbed dry before weighing. Snakes were an average of 22.1 ± 3.8% muscle by mass, but this proportion varied regionally due to uneven total segment masses and viscera masses; absolute muscle mass was highest at midbody and decreased anteriorly and posteriorly, though the highest muscle mass was only 27.0 ± 3.8% higher than the average.
Because of postmortem and preservation distortion of body shape (facilitated by the mobile ribs of snakes), average body radius was computed from SVL and mass by treating the snake as a uniform cylinder with a tissue density of 1.05 g cm3 (typical for vertebrate tissues). As the muscular lever arms for lateral flexion are unknown for any snake, we estimated the maximum lever arm as 1/2 the radius of the body; while some epaxial muscles (e.g., m. iliocostalis) may have larger lever arms, others likely have much lower lever arms (m. multifidus and semispinalis-spinalis). Similarly, muscular PSCA for the entire epaxial muscle group was computed as a cylinder based on SVL and total muscle mass; while snake epaxial musculature is highly complex, none of the muscles show strong pennation. Peak isometric muscle force was estimated based on the standard 30 N cm2 value seen in most vertebrate muscles, and divided by two to account for unilateral activation (Jayne, 1988). Although the maximal shortening speed and shape of the force-velocity relationship is unknown in snakes, we assumed that during lateral undulation, snakes would be operating near their peak isotonic power, and thus with a force of half the peak isometric muscle force; as activation/deactivation kinetics and length tension properties are also unknown in snakes, we did not attempt to account for these. Peak torque was computed as this force divided by estimated lever arms. Although many crucial properties are unknown, this value represents a charitably high estimate of peak torque; this value would be depressed by steeper force-velocity curves, departure from the plateau of the length-tension curve, incomplete activation during the work loop, or lower muscle lever arms, while the higher muscle mass at midbody and slightly larger vertebrae would increase the peak torque.
Data availability
Tracked snake midlines and granular drag data have been deposited to SMARTech online repository. Codes are available on GitHub at https://github.com/PerrinESchiebel/SnakeTrackingCodes (copy archived at https://github.com/elifesciences-publications/SnakeTrackingCodes).
-
SMARTechMitigating memory effects during undulatory locomotion on hysteretic materials dataset.
References
-
Slow drag in a granular mediumPhysical Review Letters 82:205–208.https://doi.org/10.1103/PhysRevLett.82.205
-
Granular drag on a discrete object: shape effects on jammingPhysical Review E 64:061303.https://doi.org/10.1103/PhysRevE.64.061303
-
Intelligence by mechanicsPhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365:199–220.https://doi.org/10.1098/rsta.2006.1911
-
Depth-dependent resistance of granular media to vertical penetrationPhysical Review Letters 111:168002.https://doi.org/10.1103/PhysRevLett.111.168002
-
Drag induced lift in granular mediaPhysical Review Letters 106:028001.https://doi.org/10.1103/PhysRevLett.106.028001
-
Force and flow at the onset of drag in plowed granular mediaPhysical Review E 89:042202.https://doi.org/10.1103/PhysRevE.89.042202
-
Semi-infinite plates dragged through granular bedsJournal of Statistical Mechanics: Theory and Experiment 2012:P07013.https://doi.org/10.1088/1742-5468/2012/07/P07013
-
BookBiologically Inspired Robots: Serpentile Locomotors and ManipulatorsOxford University Press.
-
What are the relations between mechanics, gait parameters, and energetics in terrestrial locomotion?Journal of Experimental Zoology Part A: Comparative Experimental Biology 305A:912–922.https://doi.org/10.1002/jez.a.335
-
The mechanics of slithering locomotionPNAS 106:10081–10085.https://doi.org/10.1073/pnas.0812533106
-
Nonlocal constitutive relation for steady granular flowPhysical Review Letters 108:178301.https://doi.org/10.1103/PhysRevLett.108.178301
-
Memory formation in matterReviews of Modern Physics 91:035002.https://doi.org/10.1103/RevModPhys.91.035002
-
Templates and anchors: neuromechanical hypotheses of legged locomotion on landJ Exp. Bio 202:3325–3332.
-
BookPhysics of Continuous Matter: Exotic and Everyday Phenomena in the Macroscopic WorldCRC press.
-
A review of fish swimming mechanics and behaviour in altered flowsPhilosophical Transactions of the Royal Society B: Biological Sciences 362:1973–1993.https://doi.org/10.1098/rstb.2007.2082
-
A review on modelling, implementation, and control of snake robotsRobotics and Autonomous Systems 60:29–40.https://doi.org/10.1016/j.robot.2011.08.010
-
ConferenceMars exploration rover mobility assembly design, test and performanceIEEE International Conference on Systems, Man and Cybernetics 2005. pp. 450–455.https://doi.org/10.1109/ICSMC.2005.1571187
-
Undulatory swimming in sand: experimental and simulation studies of a robotic sandfishThe International Journal of Robotics Research 30:793–805.https://doi.org/10.1177/0278364911402406
-
Mechanical models of sandfish locomotion reveal principles of high performance subsurface sand-swimmingJournal of the Royal Society Interface 8:1332–1345.https://doi.org/10.1098/rsif.2010.0678
-
Friction enhancement in concertina locomotion of snakesJournal of the Royal Society Interface 9:3067–3080.https://doi.org/10.1098/rsif.2012.0132
-
Flipper-driven terrestrial locomotion of a sea turtle-inspired robotBioinspiration & Biomimetics 8:026007.https://doi.org/10.1088/1748-3182/8/2/026007
-
Springer Handbook of Robotics1151–1173, Search and rescue robotics, Springer Handbook of Robotics, Springer.
-
Lattices of hydrodynamically interacting flapping swimmersPhysical Review X 9:041024.https://doi.org/10.1103/PhysRevX.9.041024
-
RHex: a simple and highly mobile hexapod robotThe International Journal of Robotics Research 20:616–631.https://doi.org/10.1177/02783640122067570
-
ConferenceReview of fish swimming modes for aquatic locomotionIEEE Journal of Oceanic Engineering.https://doi.org/10.1109/48.757275
-
Locomotor benefits of being a slender and slick sand swimmerJournal of Experimental Biology 218:440–450.https://doi.org/10.1242/jeb.108357
-
Undulatory and oscillatory swimmingJournal of Fluid Mechanics 874:284.https://doi.org/10.1017/jfm.2019.284
-
Neuromechanical response of musculo-skeletal structures in cockroaches during rapid running on rough terrainJournal of Experimental Biology 211:433–446.https://doi.org/10.1242/jeb.012385
-
Dimensionality and dynamics in the behavior of C. elegansPLOS Computational Biology 4:e1000028.https://doi.org/10.1371/journal.pcbi.1000028
-
A perspective on medical roboticsProceedings of the IEEE 94:1652–1664.https://doi.org/10.1109/JPROC.2006.880669
-
PIVlab – Towards User-friendly, Affordable and Accurate Digital Particle Image Velocimetry in MATLABJournal of Open Research Software 2:53315003.https://doi.org/10.5334/jors.bl
-
The hydrodynamics of eel swimming II. effect of swimming speedJournal of Experimental Biology 207:3265–3279.https://doi.org/10.1242/jeb.01139
-
The effectiveness of resistive force theory in granular locomotionPhysics of Fluids 26:101308.https://doi.org/10.1063/1.4898629
Article and author information
Author details
Funding
National Science Foundation (PHY1205878)
- Perrin E Schiebel
- Henry C Astley
- Jennifer M Rieser
- Alex M Hubbard
- Daniel I Goldman
Army Research Office (W911NF-11-1-0514)
- Jennifer M Rieser
- Perrin E Schiebel
- Henry C Astley
American Society for Engineering Education (32 CFR 168a)
- Perrin E Schiebel
Army Research Office (W911NF-18-1-0120)
- Perrin E Schiebel
- Henry C Astley
- Jennifer M Rieser
- Christian Hubicki
- Alex M Hubbard
- Kelimar Diaz
- Daniel I Goldman
Simons Foundation
- Perrin E Schiebel
- Kelimar Diaz
- Daniel I Goldman
National Science Foundation (PHY-1150760)
- Perrin E Schiebel
- Jennifer M Rieser
- Henry C Astley
- Alex M Hubbard
- Daniel I Goldman
National Science Foundation (CMMI-1361778)
- Perrin E Schiebel
- Jennifer M Rieser
- Henry C Astley
- Daniel I Goldman
Army Research Office (W911NF-18-1-0118)
- Ken Kamrin
- Shashank Agarwal
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Kevin and April Young for collecting the C. occipitalis and Mark Lowder for his assistance in creating the non-sand-specialist tracking code. This work was supported by NSF PoLS PHY-1205878, PHY-1150760, and CMMI-1361778. ARO W911NF-11-1-0514, ARO W911NF-18-1-0120, ARO W911NF-18-1-0118 (KK and SK), NDSEG 32 CFR 168a (PES), and the Simons Southeast Center for Mathematics and Biology.
Ethics
Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All experimental procedures using the non-sand-specialists were conducted in accordance with the Georgia Institute of Technology institutional animal care and use committee (IACUC) protocol A14058. All procedures involving these zoo animals were reviewed and approved by the Zoo Atlanta Scientific Research Committee. All experimental procedures using C. occipitalis were conducted in accordance with the Georgia Institute of Technology IACUC protocols A14066 and A14067. All experiments were non-invasive and every effort was made to minimize distress.
Copyright
© 2020, Schiebel et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 832
- views
-
- 114
- downloads
-
- 30
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Physics of Living Systems
Many proteins have been recently shown to undergo a process of phase separation that leads to the formation of biomolecular condensates. Intriguingly, it has been observed that some of these proteins form dense droplets of sizeable dimensions already below the critical concentration, which is the concentration at which phase separation occurs. To understand this phenomenon, which is not readily compatible with classical nucleation theory, we investigated the properties of the droplet size distributions as a function of protein concentration. We found that these distributions can be described by a scale-invariant log-normal function with an average that increases progressively as the concentration approaches the critical concentration from below. The results of this scaling analysis suggest the existence of a universal behaviour independent of the sequences and structures of the proteins undergoing phase separation. While we refrain from proposing a theoretical model here, we suggest that any model of protein phase separation should predict the scaling exponents that we reported here from the fitting of experimental measurements of droplet size distributions. Furthermore, based on these observations, we show that it is possible to use the scale invariance to estimate the critical concentration for protein phase separation.
-
- Computational and Systems Biology
- Physics of Living Systems
Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.