Mitigating memory effects during undulatory locomotion on hysteretic materials

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.

The generalist pygmy rattlesnake Sistrurus miliarius, which lives throughout the southeastern US, attempting to move on natural sand collected from Yuma, Arizona, USA. The animal has completed several undulations, sweeping GM lateral to the midline of the body. This snake failed to progress further than pictured. (b) The generalist eastern indigo snake Drymarchon couperi lives in the southeastern US, although its range is smaller than S. miliarus, primarily inhabiting Florida. On the same GM as (a) (c) The sand-specialist shovel-nosed snake Chionactis occipitalis in the lab on 300 µm glass particles. (d-i) Digitized midlines of animals. Color indicates time from beginning to end of the trial. All scaled to the 10 cm scale bar shown. (d-g) are on the Yuma sand and (h,i) are on glass particles. (d) Crotalus tigris. This snake superficially shares some of its natural range with that of C. occipitalis, but is closely associated with rocky boulder substrates rather than sand. The animal was unable to progress on the GM. Total length of the trial ttot =25.7 s, time between plotted midlines δt = 100 ms (e) Crotalus lepidus, a generalist in the southwestern US and the central part of Mexico. ttot =9.4 s, δt = 33 ms (f) Dasypeltis scabra inhabits a wide range of habitats in Africa. ttot =1.6 s, δt = 33 ms (g) Nerodia sipedon, a water snake inhabiting most of the Eastern US extending into Canada. ttot =6.0 s, δt = 33 ms (h) C. occipitalis 130 and (i) 128. These trials represent the intra-individual variation in kinematics. Some of the animals moved approximately "in a tube" where all segments of the body followed in the path of their rostral neighbors (e.g. (h), ttot =1.25 s, δt = 12 ms) while others used this strategy only on the anterior portion of the body, appearing to drag the posterior segments in a more or less straight line behind themselves (e.g. (i), ttot =1.01 s, δt = 12 ms). See [21] for tables of length, width, and mass of all animals used. . θm was comparable to values measured from images of tracks taken in the field [21]. (inset) shows a close-up of the data with mean and range measured in each trial plotted individually and color consistent with the main plot (Linear fit with 95% confidence interval to inset data: slope -0.016 (-0.020,-0.011) intercept 2.66 (2.42, 2.90), R-square=0.6). (c) Axes are as in (b), black cross represents the average C. occipitalis measurements. Colored crosses are the mean and range of values measured on a per-trial basis using the various non-specialist species, indicated by color. N=22, n=38.

121
The forces acting on animals moving submerged within GM were elucidated by subsurface drag measurements [32]. 122 These experiments revealed that during subsurface swimming in GM the material acts like a frictional fluid; propulsive 123 forces arise from the resistance of the GM to motion of the body segments as the grains interact with each other 124 via normal and frictional contacts [15,32]. A similar understanding of the character of movement at the surface is 125 unknown. 126 We observed the formation of granular piles as the snakes self-deformed on the granular surface ( Fig. 3(a), [21]). 127 This suggested that, consistent with subsurface swimming and surface walking, slithering animals propel themselves self-deformation pattern. 133 We empirically measured the granular stress on a simple model for a snake body segment-an aluminum plate 134 commanded to move at a slip angle β d between the drag velocity unit vectorv d and the plate face tangent in the 135 horizontal plane ( Fig. 3(a,b)). We dragged the plate for 20 cm at a constant depth, z, measured from the intruder's 136 bottom edge ( Fig. 3(b)) and measured stress normal, σ n , and tangent, σ t , to the plate face (Fig. 3(c), Appendix E).

137
A fluidizing bed containing the same 297 ± 40 µm glass particles used in the C. occipitalis experiments prepared the 138 material to an initially loose-packed state.

139
Unlike subsurface drag stresses, which developed almost instantaneously to the steady state [32], at the surface stress 140 monotonically increased over several centimeters before saturating (Fig. 3(d)). This is due to the free surface flow of 141 the GM; a pile of sand above the surface, like those created by the snakes, appeared at the leading face of the intruder 142 at the onset of drag and increased in volume until reaching a balance between the new grains being encountered and . β d = 90 • for all trials such that the total stress was equal to σn. Circle markers are mean and error is std. of three trials. Gray curve is model prediction for z = 8 mm and light blue curve for z = 12 mm. (inset) Normal force versus v d measured in experiment (circle markers) and as predicted by the model (solid curves). Color is consistent with main plot. Each dot is the average force measured in one trial, all trials shown. (b) Stress normal to the plate face as a function of z at v d = 10 mm s −1 .
A. Velocity-dependent stress captured by grain inertia model.

148
The snake speeds were variable (segment speeds from ≈ 35 to 95 cm s −1 [21]) and the intrusion depth of the snakes' 149 trunk into the GM ranged from 0 (no intrusion, occurring at the apexes of the wave which the snake lifted off of 150 the surface) to ≈ 5 mm [21]. Previous studies indicated that granular drag stress depends on both intruder speed 151 [36] and depth [26,37]. Commensurate with those studies, we fixed β d = 90 deg and z = 8 mm and found normal 152 stress quadratically increased as v d increased from 1 mm s −1 to 750 mm s −1 (the limit of the robot arm capability, 153 Fig. 4 (a)). Similarly, for β d = 90 deg and v d = 10 mm s −1 , normal stress linearly increased with z from 4 mm, the 154 shallowest depth where force could be resolved, to 40 mm, where the plate was fully submerged with the top edge 155 10 mm below the surface ( Fig. 4(b)).

156
For snakes moving slowly on non-deformable surfaces, body inertia was small compared to the frictional forces 157 between the belly and the surface [33]. While C. occipitalis was moving quickly, we observed, in line with earlier  The animal behavior indicated that friction between the body and the GM and dissipative interactions within the 162 GM had greater impact on the motion than inertia of the body. We thus assumed a resistive-force-dominated system 163 in which stress arose from independent static and dynamic terms. 164 We modelled stress on the plate as a static stress σ o (v d → 0), assumed to arise from Coulomb-wedge type forces 165 [35], plus a material inertia term αρv 2 d (GM density ρ = φρ glass = 0.68 × 2, 500 kg m −3 ). We estimated a packing Since σt is small, values with noise included can approach zero. As we were interested in the force evolution over several cm (for the plotted drag speed this corresponds to signals < 1 Hz), we removed fluctuations above 5 Hz using a low-pass butterworth filter. Curves from three trials are shown for each of four different, fixed β d = 15 (darkest, bottom curve), 30, 45, and 60 deg (lightest, top curve) as labeled on the right of the plot. The slopes of a linear regression fit to the average of three trials as a function of drag distance were −1.8 × 10 −5 , −2.9 × 10 −5 , −5.0 × 10 −5 , and −6.4 × 10 −5 cm −1 for β d = 15, 30, 45, and 60 • , respectively (R-square=0.6, 0.8, 0.8, and 0.7). (b) Average stress as a function of changing drag angle. Each measurement was calculated by averaging the raw stress data from 10 to 20 cm drag distance. Each data point is the mean and standard deviation of three trials; error bars are smaller than the markers. Solid lines are the fit functions used in the RFT calculations.
fraction φ = 0.58 induced by motion of the intruder with a scalar α which was expected to be of O(1). We estimated 167 σ o using the experimental data by subtracting ρv 2 d from σ n measured at the three lowest speeds and averaging the of drag distance were small relative to the effect of changing β d . Thus, we averaged σ n and σ t to characterize the 189 relationship between plate orientation and stress normal and tangential to its face ( Fig. 5(b)). 190 We found that normal stress increased monotonically with β d while tangential stress was approximately constant, 191 gradually falling to zero as β d went to 90 • (Fig. 5(b)). The ratio σn /σt as a function of β d for both varying speeds and 192 depths collapsed to a characteristic anisotropy curve ( Fig. 6(a)). The same drag anisotropy curves were measured in 193 poppy seeds and oolite sand and were similarly independent of depth and speed [24]. The comprehensive appearance 194 of this curve suggests it may be a more general feature of dissipative, deformable materials.

195
Because tangential stress magnitudes were relatively small and normal stresses rose sharply with β d at small angles 196 ( Fig. 5(b)), thrust and drag forces on the plate were equal ( σn /σt = 1) at smaller β d than in viscous fluid or subsurface 197 in GM ( Fig. 6(a)). Reflecting the small angles where the anisotropy was unity, β s measured on the snake were small 198 ( Fig. 6(a), gray histogram).

199
Plotting anisotropy at constant β d as speed and depth varied further illustrated the relatively weak dependence 200 of σn /σt on these variables ( Fig. 6(c,d)). Especially at small β d , anisotropy did not depend on v d . σn /σt was more 201 dependent on z at shallow depths, decreasing twofold as z increased from 4 to 8 mm (Fig. 6(c)). inexpensive. RFT successfully models a number of systems for which these assumptions are valid, although its 208 effectiveness in a system exhibiting hysteresis was unknown [15].

209
The negligible body inertia and local, dissipative nature of the forces indicated using resistive force theory (RFT,

210
[40]) and assuming quasi-static locomotion, such that forces were balanced at each moment, was a good candidate 211 for probing the connection between waveform and performance of the surface slithering snakes. The insensitivity of 212 anisotropy on speed suggested that the same speed-independent granular RFT used to predict subsurface performance 213 in GM [15, 32] could be applied to this system.

214
The empirical relationships between body segment motion, characterized by β d , and stress measured in the drag 215 experiments ( Fig. 5(b)) provided the necessary information to estimate the forces acting on different waveforms using 216 RFT.

217
We first focused on calculation of sand-specialist performance. RFT, using the force relations for v d = 10 mm s −1 218 and z = 8 mm, accurately predicted the average v CoM and the relationship between v CoM and v seg of C. occipitalis 219 for the average snake kinematics and morphology [21]. We also performed test RFT calculations using the drag 220 stresses measured at each depth and found that the effect of depth was much less than that of shape and less than increasing the tangential stress shifts the anisotropy curve down, moving the location of σn /σt = 1 to larger β d 234 ( Fig. 6(a)).

235
A. Trade off between internal and external factors constrains waveforms. 236 We used RFT to predict the maximum joint torque, τ m,RF T , experienced by any C. occipitalis body segment over as θ m or ξ decreased both because balancing forces required larger β d , thus greater force magnitudes ( Fig. 5(b)), and 239 these shapes "stretched out" the body, creating longer lever arms (see Fig. 7(b) and compare shapes in lower left 240 corner to upper right).

241
The snakes could minimize torque by increasing ξ and/or θ m , however, we did not observe these waveforms (Fig. 2). Similarly, snakes did not use waveforms which minimized mechanical cost of transport or distance traveled per cycle 243 [21]. We thus hypothesized there were internal factors which "penalized" the low-torque waveforms.

244
If, for simplicity, we assume the animal moved with no slip regardless of waveform such that the distance traveled 245 per cycle was equal to the wavelength, there were two (not necessarily independent) ways to increase v CoM : (1)  In limbed organisms, movement speed is related to the interplay between gait parameters like stride length and 249 frequency and physiological concerns like energetic cost (e.g. [42,43]). We hypothesized that, because the snake 250 performed an escape response in our experiments, their objective was to maximize speed.

251
For a desired v CoM , we explored how quickly a "muscle" segment of fixed nominal length relative to total body 252 length would have to shorten as a function of ξ and θ m . The shortening speed was a function of both how quickly 253 curvature changed along the body for a given waveshape and the frequency needed to achieve the target v CoM for that 254 shape's stride length. We modeled a body segment as a bending beam with length along the spine δs and arclength 255 along the inside of the curve δs (see diagram Fig. 7(b)). We approximated the amount of shortening the muscles 256 must produce as the difference, ∆s = δs − δs , at the point of maximum bending. Using this model we estimated 257 v shorten = δs / 1 4 T , the speed the inside of the bend must change length to go from δs = δs at the inflection points of the 258 wave to the maximum δs occurring at the peaks given undulation period T set by the frequency (see Appendix G).   C. occipitalis performed an escape response in our trials. Thus, we hypothesized that the specific waveform used 272 by the animals was that which endowed maximum speed across the substrate. We endeavored to include the trade off 273 between decreasing torque and increasing actuation speed needed to maintain v CoM (now using RFT to estimate the 274 actual stride length given slipping) as θ m and ξ increased. To do so, we used RFT to calculate the joint-level power 275 (the rate of change of the joint angle dotted with the torque experienced by that joint, Appendix F) of each segment 13 over a cycle for each shape in the θ m , ξ space. The power-limited velocity, v pl , was the CoM speed for which the peak 277 power generated by any joint over the cycle was equal to a constant, peak available power.

278
In accord with the tradeoff between internal actuation demands and external torques, v pl was maximized in the 279 center of the θ m , ξ space inhabited by C. occipitalis' waveform ( Fig. 8(a)). Reflecting the oscillations in τ , this metric 280 had maxima near integer values of ξ.

281
Power-limited velocity was maximized at the number of waves (Fig. 8(b)) as well as at the attack angle (Fig. 8(c)) 282 used by C. occipitalis. We estimated the peak torque output of C. occipitalis muscles using dissection of museum 283 specimens (Appendix I). The attack angle used by the animals was near the point where τ m,RF T increased above the 284 estimated muscle capability. We report both maximal torque output (Fig. 8(c) upper horizontal bar) as well as output 285 scaled for the average contraction velocity estimated as the body-wall shortening speed (Fig. 8(c) lower horizontal 286 bar). It may be that the individual variation in θ m reflects differences in peak muscle power capabilities.

287
C. Stout snakes must use larger attack angles to succeed.

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

291
A striking difference between the various species was their aspect ratio, L /w, the ratio of the total length to the 292 diameter of the body at the widest point. A slender snake like Fig. 1(b,c) will thus have a higher aspect ratio than 293 a stout one, e.g. Fig. 1(a). We previously discovered that C. occipitalis' high-aspect-ratio allowed it to move more 294 effectively when swimming subsurface in GM than a low-aspect-ratio sand-swimming lizard [5]. 295 We characterized the performance of the snakes using the unsigned average slip angle, β s = acos(|t ·v|) (  We found that animals with β s > 30 degrees were those that were ineffective ( Fig. 9(a), red crosses). Some snakes 301 with slip close to 30 deg were able to make consistent forward progress (e.g. Fig. 1(e)) while the tracks of those 302 animals with β s > 30 deg either did not progress at all (e.g. Fig. 1(a,d)), or would have only occasional spurts of 303 forward motion linked by extended periods of undulating in place. 304 We observed a transition in the slip versus aspect ratio plot at an aspect ratio of about 26 ( Fig. 9(a), vertical 305 dashed line). At larger L /w the snakes were, with the exclusion of one exception, moving effectively with an average 306 β s = 10.3 ± 3.8 deg ( Fig. 9(a), blue crosses). Below L /w = 26 performance was highly variable. Four of the five species 307 which failed to progress across the level sand had L /w < 26 ( Fig. 9(a), red crosses). The other snakes of L /w < 26 308 moved with higher slip on average than their high-aspect ratio counterparts, β s = 20.3 ± 5.8 deg ( Fig. 9(a), gray 309 crosses). This trend was recapitulated in RFT calculation of slip as a function of aspect ratio for a snake of fixed 310 waveform and mass ( Fig. 9(a), [21]).

311
RFT predicted that as attack angle decreased, the slip of a snake with fixed body morphology and wavenumber 312 would increase ( Fig. 9(b), black line). We examined slip as θ m changed for those snakes with an aspect ratio less than 26 (and the snake that failed with L /w=30). Those which failed to progress ( Fig. 9(b) red crosses) were at lower attack 314 angles than those which succeeded ( Fig. 9( had L /w of greater than 26 generally used higher attack angles as well ( Fig. 9(b, inset) blue crosses).

317
RFT was reasonably accurate in predicting the performance of snakes with β s < 30 deg, given the average θ m , ξ, 318 length, and mass of the individual ( Fig. 9(a), black and green squares). However, RFT underestimated the slip of 319 those snakes which failed ( Fig. 9(a), red squares). 320 We observed that the snakes created permanent disturbances in the surface of the GM. While those which moved 321 with low slip created observable piles of sand at the posterior-facing side/s of the body ( Fig. 1(b,c)), those which failed 322 appeared to dig themselves into a channel, primarily depositing GM laterally to the long-axis of the body ( Fig. 1(a)).

ζ= δθ δs
Markers 10. Systematic exploration of memory effects using a robophysical model. (a) 10 joint, 11 segment robot. The shape was a serpenoid curve, θ(n, t) = θmsin(2π(ξ n N + f t)) for joint n out of N total joints, controlled by commanding the servo motor positions, ζ, to vary sinusoidally in time with a fixed amplitude. Each motor was offset from its anterior neighbor by a constant phase set by ξ. Length=72 cm and L /w = 11.3. (b,c) Stills of the robot taken after it has completed 2.75 undulation cycles. θm = 52 deg in (b), θm = 105 deg in (c), and ξ = 1 in both cases. Robot is artificially colored yellow to distinguish the body from the tail cord. Note the large appearance of the tail cord in c is due to perspective as the researcher was standing and loosely holding the cord above the substrate.

324
The waveform and performance of the snakes was variable, and this variability was reflected in the tracks left by 325 the animals. As the material did not re-flow around the body after being disturbed, we hypothesized that the manner 326 in which different waveforms remodeled the substrate was important in determining performance. We systematically 327 explored the impact of waveform on performance using a robophysical model, a 10 joint robot on the surface of poppy 328 seeds (Fig. 10).
We commanded the motors at each joint of the robot to execute a serpenoid curve (θ(s, t) = θ m sin(2π ξ /L(s + v seg t))

331
[2]). We varied ξ and θ m and recorded kinematics as the robot performed three undulation cycles on the surface of 332 poppy seeds (Appendix H). The number of joints on the robot was limited by motor strength; as a result the robot 333 could only achieve ξ between 1 and 1.4 as lower ξ required too much torque and higher were not well-resolved by the 334 number of joints. The power per unit volume of the motors also limited the aspect ratio; L /w = 11.3 is lower than 335 that of the stoutest snakes tested. 336 We characterized performance by measuring slip as in the snakes (Fig. 11(a)). We also measured the CoM displace-337 ment in a single undulation cycle normalized by the total length (BLC, Fig. 11(b)). This variable provided intuition 338 for how effectively the robot was progressing across the substrate.

339
Performance of the robot was a function of the waveform parameters. This was intuitive given our RFT predictions 340 for the snake, which depended on both ξ and θ m (e.g. Fig. 8), and was commensurate with our findings in the biological 341 experiments ( Fig. 9(b)). Over the range of ξ available to the robot performance was more strongly dependent on 342 attack angle than wavenumber, therefore we focused our attention on variations in θ m .

343
B. RFT prediction of robot performance is inaccurate. 344 We used granular stress relations measured in poppy seeds using a flat plate of the same plastic used to print the 345 robot body [24] in a resistive force theory calculation predicting the robot performance as a function of θ m for fixed 346 ξ = 1 (Fig. 11(a,b), black curves).

347
In previous studies, RFT accurately predicted locomotor performance of a number of systems [15] and in our study 348 RFT was accurate for C. occipitalis as well as the successful non-specialist snakes. However, as in the high-slip 349 snakes, RFT under-predicted slip at low attack angles (Fig. 11(a)). Despite the substantial predicted slip, RFT 350 indicated these waveforms would still make forward progress ( Fig. 11(b,e)). The robot, however, did not displace at 351 the smallest attack angles tested ( Fig. 11(b,c)). Conversely, at high θ m RFT under-predicted displacement of the 352 robot ( Fig. 11(b)).

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

358
Material hysteresis, which was not accounted for in our RFT calculation, appeared to play an important role.  The material flow was a function of the waveshape. We characterized the overall structure of this flow by measuring 361 ψ v , the angle between the velocity vector of the grains at a given point in time and space and the overall direction of 362 motion of the robot (Fig. 12(a)). We estimated grain velocities using Particle Image Velocimetry (PIV, Fig. 12(a)), Appendix H) and calculated the probability density to measure a given ψ v over a run for ξ = 1 and each of the three 364 θ m tested ( Fig. 12(b)).

365
At high attack angles significant amounts of the material were deposited by the posterior-facing body segments, 366 similarly to what we observed in the successful snakes (Fig. 10(c), [21]). Consistent with this observation, we measured 367 peaks in the probability density of ψ v near ±180 deg, corresponding to grains which are flowing opposite the average 368 direction of motion of the robot (Fig. 12(b), orange curve). When the granular piles moving with the body at 369 later points in time re-enountered these piles the additional stress resisted backward slipping of the robot, increasing 370 displacement above that predicted for a frictional fluid (Fig. 11(b)).

371
The GM being pushed by low θ m waveforms had a significant velocity component lateral to the direction of motion; 372 as θ m decreased, the peaks in the ψ v probability density shifted toward ±90 deg (grain flow perpendicular to the 373 average direction of motion, Fig. 12(b), green and blue curves). These piles were thus deposited to the side of the 374 robot and with each undulation more material was added to these "sidewalls" which can be seen in Fig. 10(b). As the GM being pushed by subsequent undulations encountered these pre-existing walls the robot was apparently unable 376 to overcome the inertia of the grains in the previously created piles, instead depositing more material in the walls 377 without changing their location. The result was that the robot swept out a trough in the GM, completely stalling 378 forward progress and sometimes even moving backward as the trunk "rolled" along the trough walls. These failures 379 were similar in appearance to those observed in the living snakes ( Fig. 1(a,d), [21]).

380
In limbed systems, RFT was often useful in predicting the performance of the first gait cycle, before material 381 re-interaction could occur [23,24]. In the robot, however, RFT was incorrect even for the first cycle (Fig. 12(c)). Enhancement or degradation of the performance was occurring within the first gait cycle as portions of the body 383 contacted GM disturbed by other segments at a previous time.

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

388
This raises an interesting subtlety: the tail is always moving through material disturbed by the head, and placing a 389 successful robot waveform back in its own tracks does not decrease performance [21]. It is apparently the cumulative 390 effect of material deposition coupled with the waveform-dependent pattern of granular reaction forces that conspires 391 to change performance.

392
Despite the complicated interplay between the robot and the substrate we found the resulting performance was 393 repeatable. For example, for all three attack angles tested at ξ = 1, each of the three trials for a given waveform 394 pushed material in a similar way and displaced a similar amount in each of the three undulations ( Fig. 12(b,c)). This 395 indicated that, while the material remodeling was complex, it was deterministic.

397
Although RFT did not consistently predict the robot's slip, it did provide a heuristic for success and give some 398 intuition as to why a locomotor with the morphology of the robot is challenged by GM. The line in the θ m , ξ space 399 which RFT predicted to be at β s ≈ 30 deg separated the successful and unsuccessful robot waveforms (Fig. 13(a)). 400 We found there were a limited range of waveforms predicted to both have β s ≤ 30 deg and be accessible to the robot 401 given its aspect ratio and joint resolution (Fig. 13(a)). In contrast, C. occipitalis had a large space of parameters for 402 which slip was small (Fig. 13(b)). This suggests that, in contrast to the robot which was sensitive to waveform, the 403 snake was robust to errors in the motor pattern.

404
The snake also had the benefit of greater resolution than the robot; the animal has ≈ 120 vertebrae (note that in 405 snakes there is variation in vertebral number between species and between individuals of the same species) while the 406 robot has only 10 joints. This flexibility allowed the snake to operate well away from low-performing shapes (low 407 ξ and θ m ). Increasing resolution would allow the robot to use more waves on the body, making a greater range of 408 low-slip waveforms available (those in the upper right-hand corner of Fig. 13(a)). However, the robot is limited by   [29] and our robot was executing predetermined waveforms. The ability to sense and respond to induced changes in 431 the substrate is apparently unnecessary given the initial selection of an appropriate pattern of self-deformation.

432
The interrelationship between morphology, waveform, terrain remodeling, and performance rationalizes why it has tail tap would elicit an escape response. If an individual did not respond to this stimulus they were returned to the 482 holding container. If an individual failed to perform for three trials in a row they would be retired from the days 483 studies. Snakes were tested at most every other day with a maximum of two successful trials collected per day. A 484 run was included if the snake performed at least four complete undulations moving along a straight trajectory at 485 apparently constant speed. 486 We tracked the black bands on the snake to determine the body midline positions using the method described in [5].

520
Given uncertainty in the surface we could not determine whether the animal was in contact with the surface or not.

521
However, we assumed that lifted segments are applying less load to the material and the track depth measurements 522 suggest that the apexes are not in contact.

523
Appendix D: Laser line measurements of the free surface 524 We characterized the tracks created by C. occipitalis using a laser-line apparatus. An off-axis camera (Logitech 525 C920 HD Pro Webcam) captured images of the laser line. Both the laser and camera were affixed to a linear bearing 526 and moved using a linear actuator (Firgelli). We used MATLAB to move the actuator and collect an image every We observed that as a body segment pushed laterally against the sand it built a pile similarly to the pile that 542 evolved during drag of the plate. The rostral-facing side of the body, opposite the pile, appeared to be in contact with 543 very little material. Therefore, we chose to model the snake body as a flat plate, representing the caudal-facing side 544 of the body, with an added term to account for drag on the ventral surface. This is not be an unreasonable picture 545 of this species, as the ventral surface is known to be flat or even slightly concave [31]. As the snake held the head 546 slightly raised from the surface we did not include a term accounting for the drag acting against the head as it is 547 pushed through the material.

548
Previous research [33] found that the ventral scutes are anisotropic such that gliding directly forward produces less 549 drag than sliding laterally and both forward and lateral motion results in less frictional drag than moving directly 550 backwards. The coefficient of static friction of C. occipitalis scales is low (0.109±0.016 ventral forward and 0.137±0.018 backward [5]). The lateral coefficient of friction for C. occipitalis scales is not known, however, given [33] results we 552 assumed that it was bounded between 0.10 and 0.14 such that the difference between frictional forces acting tangential 553 and normal to segments was negligible compared to the larger plate-drag forces. Therefore, we chose not to include the 554 frictional anisotropy in our RFT calculation. We previously found the coefficient of static friction between Aluminum 555 and the glass particles is approximately 0.2 [32]. Therefore, in predicting the performance of C. occipitaliswe scaled 556 the σ t function measured using the Aluminum plate by 1 /2.

557
In modeling the various snake species we used frictions 0. We modeled ventral drag on each segment as acting opposite the velocity with magnitude µmg/n segs (µ = 0.1 was 562 the coefficient of static friction for the snake's ventral scutes on the glass particles [5], m was the average snake mass, 563 g = 9.81 m s 2 was the gravitational constant, and n segs the number of segments used in the calculation).

564
Given the low-friction scales of C. occipitalis, including lifting in the calculations did not substantially improve 565 predicted performance [21]. We speculate that these segments may be lifted as a side-effect of the muscle activation 566 responsible for generating the horizontal waveform, that is, the morphology of the trunk means that vertical lifting is 567 a consequence of generating the curvature in the horizontal waves. It may also be that the lifting is intentional and a 568 buffer against deleterious motor program mistakes or changes in the environment, or the small benefit of lifting these 569 segments is greater than the energetic cost. We did not re-distribute ventral drag forces to account for lifting. 570 We approximated the kinematics of the snake using 100 segments at 70 points in time divided evenly over one full 571 undulation cycle. These values are well past those needed for convergence of the calculation. We calculated each 572 segment's orientation and velocity using θ(s s , t j ) = θ m sin( 2πξ /L(s i + v seg t j )) for each combination of parameters θ m 573 and ξ. Summing over the body yielded the total force acting on the CoM in x and y and torque about the CoM.

574
MATLABs lsqnonlin function was used to find the x and y components of v CoM and the angular velocity about the 575 CoM which resulted in zero net stress. 576 We calculated torque acting on the CoM using i r i × F i where r i is the location of each segment with respect to the CoM and F i are the RFT calculated forces acting on that segment. We calculated torque acting at each joint as in [6]. For joint number n the internal torque is τ n,int = n i=1 ( r i − r n ) × F i . To predict τ m,RF T in a cycle we find this value for each of the segments at all times. τ m,RF T is the maximum value in this 100 × 70 matrix. Joint-level power for each joint at each time was calculated using the rate of change of the joint angle, ζ n = dθn ds , and power was thus We kinematically model a snake body segment as a beam with midline length δs and length of the inside of the curve δs . For a given body radius, r, and radius of curvature, R, δs = δs R−r R . The speed the inside of the segment must change length from the unbent length of δs to δs is v shorten = δs−δs δt . We can write δs − δs as δs(1 − R−r R ). Using the relation R = κ −1 we get v shorten = δsrκ. Next, using κ = dθ ds for a serpenoid curve θ(s, t) = θ m sin(2π( ξ L s + f t)), we find the maximum curvature in terms  Because of postmortem and preservation distortion of body shape (facilitated by the mobile ribs of snakes), average 597 body radius was computed from SVL and mass by treating the snake as a uniform cylinder with a tissue density of 598 1.05 g cm 3 (typical for vertebrate tissues). As the muscular lever arms for lateral flexion are unknown for any snake, 599 we estimated the maximum lever arm as 1 /2 the radius of the body; while some epaxial muscles may have larger lever 600 arms (m. iliocostalis), others likely have much lower lever arms (m. multifidus and semispinalis-spinalis). Similarly, 601 muscular PSCA for the entire epaxial muscle group was computed as a cylinder based on SVL and total muscle 602 mass; while snake epaxial musculature is highly complex, none of the muscles show strong pennation. Peak isometric 603 muscle force was estimated based on the standard 30 N cm 2 value seen in most vertebrate muscles, and divided by 604 two to account for unilateral activation [52]. Although the maximal shortening speed and shape of the force-velocity 605 relationship is unknown in snakes, we assumed that during lateral undulation, snakes would be operating near their 606 peak isotonic power, and thus with a force of half the peak isometric muscle force; as activation/deactivation kinetics 607 and length tension properties are also unknown in snakes, we did not attempt to account for these. Peak torque was 608 computed as this force divided by estimated lever arms. Although many crucial properties are unknown, this value 609 represents a charitably high estimate of peak torque; this value would be depressed by steeper force-velocity curves, 610 departure from the plateau of the length-tension curve, incomplete activation during the work loop, or lower muscle 611 lever arms, while the higher muscle mass at midbody and slightly larger vertebrae would increase the peak torque. waveform. There were two dominant PCs which were well-fit by sinusoids with phase difference π/2 ( Fig. S1(b)).

716
A plot of α 1 versus α 2 revealed the coefficients trace out a circle in time (Fig. S1(c)) which, when combined with the 717 sinusoidal PCs, corresponds to a traveling sinusoidal tangent angle along the body. Such a wave is called a serpenoid  object depend primarily on cross-sectional area with very little dependence on the shape of the object, much less than 731 in fluids [3]. While there is a shape-dependence on forces in the vertical direction [4], we assumed that vertical forces 732 acting on the snake were balanced and did not affect the horizontal forces determining performance. surface RFT to calculate mCoT as θ m and ξ varied (Fig. S7(a)). mCoT increased when θ m or ξ were small as to 737 achieve zero net stress these shapes required segments to be at larger values of β s therefore yielding more material in v CoM per increase of v seg . BLC decreased as the number of waves on the body increased ( Fig. S7(b)). The 746 body length was fixed, so adding waves on the body shortened "stride length" Changing ξ did not change mCoT as 747 drastically because decreasing the arc length of a single wave (by increasing ξ) decreased both power and v CoM . The 748 animals appeared to be using a waveform which increased BLC in comparison to the subsurface shape (crosses in 749 ( Fig. S7(b)), however, BLC was maximized at a lower ξ than used by the animal. 750 We did not include in our calculation the energy required to overcome internal resistance from viscous and/or elastic where N=100 is the number of body segments and w int is the width of the Aluminum plate used in drag experiments. 760 We assumed that the intrusion depth of the snake body was equal to that of the plate, z = 8 mm.     (a). (c) Probability density of the wavenumber of the horizontal wave (gray) and vertical wave (red). The ratio of the median values of the two curves is 2.1, indicating that the spatial frequency of the wave of lifting is twice that of the wave in the horizontal plane. There were some outliers due to measurement error which were ξ > 10 that have been excluded as unphysical. (d) The phase of the wave in z versus the wave in θ. The wave in z is traveling at twice the frequency of that in θ. The result is that the peaks of the wave in z align with the extrema in θ, that is, the apexes of the horizontal wave are lifted off of the substrate. N = 3 individuals (122-11 trials, 124-3 trials, 131-3 trials) Ratio of the median values of each curve is 2.1 An anisotropy factor of two corresponds to the animal, as the low-friction scales of the snake reduces the friction by a factor of two as compared to the Aluminum plate used in the drag measurements. The blue line is RFT calculation for the snake parameters without any lifting, light blue calculated for lifting of segments which are in the top 38% of curvatures, and green dashed for in the top 73%. The black line and gray bar is average and range of vCoM calculated from experiment. The RFT prediction and animal measurements come into agreement at an anisotropy factor of two, as expected. We also see that lifting does not greatly impact vCoM , nor does decreasing friction (thereby increasing the anisotropy factor) beyond the value achieved by the snake. (b) RFT calculation of vCoM as a function of θ. The different colors indicate force relations obtained by fitting data taken at the noted depths. We find that the RFT prediction was relatively insensitive to depth which we expected given the depth independence of K. Further, the peak of the curves does not depend on depth. Therefore, we use in the calculations force relations measured at a depth of 8 mm.