Locomotion (or self-propulsion) is an essential behavior in most living systems (1,2) and important for engineered devices like robots (3)(4)(5). In organisms as diverse as jumping kangaroos, swimming eels, crawling nematodes, and spiraling bacteria, self-propulsion results from cyclic changes in body and/or appendage configuration. These configuration sequences are ultimately generated by numerous interacting and coordinated components coupled to environments of varying composition. A major challenge in locomotor biology is to discover general principles that govern how organisms generate and control fast, stable, or energetically efficient locomotion. At the organismal scale, these principles have long been discussed in the physiology and motor control literature (6,7) with the term "neuromechanics" used to indicate the importance of concomitant consideration of nervous, musculoskeletal, and biomechanical systems in explaining performance.

Several approaches are used to develop neuromechanical control principles. One approach ("bottom-up") to addressing this question is to directly incorporate the nonlinearly coupled nervous/musculoskeletal systems and environments in full detail. While success has been achieved using this approach in ferreting out mechanisms of legged (8,9) and undulatory locomotion (10)(11)(12), the complexity of such models leads to challenges in discovering broad principles. Another approach ("top-down") ignores the complexity of organisms and seeks to discover broad (cross-taxa) and relatively simple patterns of dynamics. These models are often referred to as templates (8,13), defined as a behavior that "contains the smallest number of variables and parameters that exhibit a behavior of interest." This approach has the benefit of producing models which are analyzable, can be used to test lower-level mechanisms, and yield insight into features across organisms. The template approach has been useful in rationalizing locomotor performance and control across taxa in legged and undulatory systems (8,(14)(15)(16)(17)(18). In addition to descriptive power, templates can also serve a prescriptive role, generating dynamics which are targets of control for the neuromechanical system that can yield beneficial locomotor properties (e.g., speed, energetics, stability), enabling robots with

A challenge in locomotion studies is to discover how organisms control their many degrees of freedom to move in diverse environments. Insight into such control schemes can aid biology and advance robots. Here, we use a geometric framework to model biological locomotion in a "top-down" approach, allowing comparison of highly damped systems across scales, levels of complexity, and environmental interactions. From tiny nematode worms locomoting in fluids to snakes and lizards in granular media, animal kinematics, and locomotor performance coincide with predictions that nearly maximize a geometric phase in the space of body configurations. This leads to a principle for forward planar undulatory locomotion in highly damped environments as well as gives insight into turning behaviors.

performance approaching those of living systems (8,(19)(20)(21). Finally, a top-down approach can offer insights into how locomotors adapt their templates to changes in locomotor morphology [e.g., lizard limblessness (22)] and/or environments [e.g., from swimming to walking in amphibious salamanders (23)].

Despite the apparent simplicity in the template-based approach, a question arises: given all the ways an organism could organize its neuromechanical system to self-deform (24)-e.g., bounce across the ground like a pogo stick (17) or send waves down its body from head to tail-how can one determine "good" ones? That is, does a general framework for locomotion exist that could provide a priori useful template dynamics which could then be used to evaluate performance and test optimality (7) in terms of speed, energy use, or stability? Surprisingly, an answer to this question arose starting in the 1980s through the work of physicists and control theorists (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36). These researchers developed a scheme which now goes by the name "geometric mechanics" (and which we will refer to as the geometric phase approach) which first uses environmental models to link small self-deformations around each body configuration to small translations and rotations in world space. Line integrals over closed paths in a space of body configurations, which represent cyclic sequences of body/limb changes, can then lead to translation/rotation in the world. These translations or rotations can be expressed as "geometric phases"-net changes in global quantities like position that depend on the shape of a cyclical path whose local parameters return to their initial values (25,26,28,29,37). These geometric phases are independent of the specific temporal details, such as the speed at which the paths are executed. In other words, the geometric phase depends solely on the path's shape and is unaffected by the specific timing or speed of the movements. Geometric phases appear in diverse situations including the Ahranov-Bohm effect, parallel transport of vectors on curved manifolds, a swinging Foucault pendulum and polarization changes of light in coiled optical fibers. * A major advance in the possibility of using geometric phase in realistic locomotor situations occurred with the introduction of the minimal perturbation coordinate (38,39). Such coordinates mitigate issues associated with the noncommutivity of translations and rotations in the plane and allow line integrals in the configuration space to be approximated by surface integrals over certain functions. These functions are referred to as "constraint curvature functions" or "height functions." Critically, height functions replace potentially laborious calculations used in the line integral approach. For example, even in simple artificial systems (40)(41)(42) to identify parameters that result in optimal performance requires considerable computational effort, comparing movements arising from an infinite combination of shape change sequences. Height functions instead enable a comparatively simple, diagrammatic approach. Their key utility is that they simplify the inverse problem: providing ready identification of gaits that maximize performance in diverse systems. Height functions also give a geometric rationalization of the marginal benefits that result from changing/adapting self-deformation patterns, without the need for significant calculation. Because of its utility, over the last decades, researchers have developed the theory so that it is applicable to a broad range of situations and applied the scheme to optimal control of artificial devices including satellite reorientation (43), robot swimming (39), sidewinding (44) and walking (22,45) in granular and frictional environments.

A particular regime of self-propulsion which could be amenable to geometric analysis is that in which dissipative forces dominate inertia. Here, cyclic patterns of undulatory selfdeformations solely dictate performance (provided the environment is uniform like in open fluid)-unlike in inertia-dominated systems where gliding (movement without shape changes) and stored/returned elastic energy can be utilized. This is typically thought of as the world of very small scales [nicely narrated in ref. 40 and the subject of much effort devoted to locomotion (46)]. Surprisingly, in the last decade, our experimental and granular resistive force theory [RFT, first introduced for microscopic organisms (47)] modeling studies of sand-swimming organisms have revealed that the dynamics of such terrestrial macroscopic undulatory locomotors (48)(49)(50)(51) operate in a mechanically analogous regime, where rate-independent friction, as opposed to viscosity, dominates inertia. Even within these dissipative locomotor regimes, comparing undulatory locomotor performance for different dynamics is challenging due again to infinite combinations of self-deformation sequences, necessitating the use of the height function formulation of locomotor geometric phase to establish cross-system principles of locomotion.

In this paper, we demonstrate that the geometric phase approach of locomotion provides a useful way to compare living organisms with seemingly very different and differently composed (e.g., exo-vs. endoskeletal) locomotor systems across scales. We first show that diverse undulatory organisms (including microscopic nematodes and macroscopic snakes/lizards, Fig. 1A) moving forward can be described using a planar wave templatea time series of two spatial basis functions. We notice that these time series exhibit similar patterns in locomotor systems across scales, suggesting the existence of a fundamental template for undulatory systems. Moreover, the time series emerge as circular paths in the two-dimensional spatial-basis space, which we refer to as serpenoid templates. Interestingly, we note that the parameters governing the animals' chosen serpenoid templates nearly maximize geometric phase, indicating that the animals are controlling their self-deformation patterns to achieve "good" A B C D Fig. 2. Low dimensional representation of animal movement. (i) Photos of animals and (ii) snapshots of animal body configurations colored by time (over two or three cycles) for (A) the nematode worm (C. elegans) in Sbasal buffer (T cycle ≈ 1 s) (B) the sandfish lizard (S. scincus) 7.6 cm below the surface of and fully immersed in 300-μm glass particles (T cycle ≈ 0.4 s), (C) the shovel-nosed snake (Ch. occipitalis) 7.6 cm below the surface of and fully immersed in 300-μm glass particles, and (D) Ch. occipitalis moving on the surface of 300-μm glass particles (T cycle = 0.3 s). (iii) Solid lines show the two dominant relative curvature ( s ) PCA modes account for (A) 96.7%, (B) 94.7%, (C) 57.6%, (D) 90.7% of the variation in observed body configurations. Dashed lines show best fits to sin and cos functions. (iv) 2D probability density map of projections of curvatures (directions specified by arrows) onto the two PCA modes with the largest eigenvalues. Axes are identical in (iii) and (iv).

locomotor performance (a notion we will discuss in more detail). Beyond forward crawling, serpenoid templates can be modulated (e.g., adding an offset to the center of the circle) to explain turning behaviors, such as the omega turn in nematodes and the differential turn in sidewinders.

Previous studies of highly damped locomotors using planar undulation revealed relative simplicity in wave shapes (52)(53)(54). To describe the planar † shapes used during C. elegans locomotion, (52) used principal components analysis (PCA) to diagonalize the covariance matrix of local body curvatures (determined from digitized midlines of animals during movement) to identify a set of orthonormal basis functions, referred to as principal components (PCs), whose weighted superposition can be used to capture observed animal body configurations. Eigenvalues associated with each PC or "eigenworm" indicate the variance explained by each mode (and therefore the fraction of the variance explained by each mode is given by the eigenvalue divided by the sum of all eigenvalues). For steady forward crawling, ref. 52 found that two PCs were sufficient to capture most of the shape variance, and the dynamics of movement could be represented within this two-dimensional space by projecting time traces of local body curvatures onto these PCs. A depiction of a "shape space" spanned by two undulatory PCs for forward movement is shown in Fig. 1B, with an example of the dynamics of movement represented by the directed closed path within this space. Each point along this path corresponds to a body posture, and the direction indicates how postures change. Environmental reaction forces induced by these posture changes (Fig. 1C) can result † Body undulation in these animals is planar but occurs dorsoventrally.

in world-frame displacements. The body shapes associated with five points identified along the path, as well as the resulting displacements, are shown in Fig. 1D.

Here, inspired by the similarity of the wave kinematics used by the sandfish lizard, S. scincus (Fig. 2B and ref. 50), and the shovel-nosed snake, Ch. occipitalis (Fig. 2 C and D and refs. 50 and 55) in sand to those used by low Reynolds number mmscale locomotor C. elegans (Fig. 2A), we investigated whether a low-dimensional representation could capture the body postures and dynamics observed in these undulatory locomotors across scales. We collected data on C. elegans (Materials and Methods) and reanalyzed previously published data on lizards and snakes (50,55). We started with the digitized animal midlines from high-speed kinematic data (Fig. 2-ii and SI Appendix, section 3), and characterized instantaneous body configuration using the relative curvature, (s, t) s , where (s, t) is the local curvature at position s and time t (Fig. 2C-ii), s is the position along the body, and s is the arc length of one wave (SI Appendix, section 3).

(s, t) s is a nondimensional and coordinate-invariant quantity that is measured as a function of position along the body for each moment in time.

PCA was applied to the entire dataset of each species combining curvature measurements from all trials throughout all times (52,56). We find that, for the forward movement of lizards and snakes in granular media, and nematodes in fluid, two PCs capture most of the variation in the body configurations of each species (SI Appendix, section 3). This observation allows us to use the space spanned by the first two PCs as our low dimensional representation for each animal:

where v 1 (s) and v 2 (s) are principal components identified from PCA; w 1 (t) and w 2 (t) are the time series of weights associated with the corresponding principal components. We thus define the shape variable (t) = [w 1 (t), w 2 (t)].

Similar to the results from ref. 52, we find that the two dominant PCs in the lizards, snakes, and nematodes were well fit by v i (s) = sin (2 ns/L + s i ), i ∈ {1, 2} where n, s i and L are, respectively, the number of waves, spatial phase, and the body length (Fig. 2-iii for modes and fits ‡ ). Note that to enforce the orthogonality of fitted basis functions, we assume | s 1 -s 2 | = /2. We use the fitted sinusoidal basis functions to construct a two-dimensional shape space. Because the first two PCs explain comparable amounts of the variance, we assign them to the dimensions of the space in order of phase (i.e., s 1 > s 2 ), rather than the amount of variance explained. This allows the chirality of the shape space to be the same across organisms (i.e., counterclockwise/clockwise directed trajectory in the shape space leading to forward/backward displacement respectively; see Fig. 2).

The visualization of shape variables reveals that animals investigated in our manuscript used nearly circular templates (of radius

to transition through body configurations in [w 1 , w 2 ]-space as they deform (Fig. 2-iv). Notably, the circle in shape space and thus the sinusoidal variation in curvature was first studied in the context of snake locomotion and limbless robot control by Hirose in his seminal work (57). Following Hirose's terminology for such a wave, we will refer to the circular pattern in shape space as a "serpenoid" template. §

To develop the geometric phase framework to rationalize such templates as well as to discover novel behaviors, we first require a model of body-environment interactions that connects a particular change in shape to center-of-mass translation and/or rotation in the environment. Resistive force theory (RFT) modeling has had successes describing movement in fluids, for example in predicting forward swimming speeds of C. elegans (59,60); however, in other situations (e.g., waves with higher curvatures), more elaborate schemes are often required to accurately capture the dynamics (61). One of the key assumptions in RFT modeling is that environmental disturbances induced by the movement of a swimmer are sufficiently localized that forces and flow fields from neighboring body segments are completely decoupled. In the last decade, studies of animal locomotion in dry granular media have revealed that the simplest form of RFT is remarkably successful in describing movement in such environments (49)(50)(51).

Further, from the assumption of decoupled forces along a deforming body, a swimmer can be divided into many infinitesimal segments that can be treated independently. In dissipationdominated environments, the net force on a body is zero at every moment in time, giving

d F ⊥ and dF are the environmental reaction forces acting perpendicular and parallel to the surface of an infinitesimal segment of the body as it moves within the surrounding medium.

‡ For the subsurface movement of Ch. occipitalis, we attribute the points near the origin to turning behavior that is not captured by the first two modes; however, we will show that despite the more complex body postures and locomotion, we are still able to quantitatively describe forward locomotion.

§ The serpenoid dynamics differ from characterization of animal body shapes as sinusoidal amplitude displacements of away from the midline of a straight animal during a posteriorly traveling wave, previously utilized in refs. 49 and 58.

Environmental reaction forces experienced by the undulatory locomotors in this study are shown in Fig. 3.

In the fluid-swimming nematodes analyzed here, locomotion occurs at sufficiently low Reynolds number (the ratio of inertial to viscous forces is approximately 0.1) such that the assumption of zero net force (inertialess locomotion) is well justified. This has an important locomotor consequence: If a nematode stops self-deforming, its locomotory speed will decay to one-half of its steady-state speed within approximately ≈5 ms via viscous Stokes drag (see SI Appendix, section 4 and ref. 62 for details on the calculation). We will refer to this as the coasting time, coast , and can form a nondimensional parameter, the "coasting number" C = 2 coast / cycle , the ratio of this time to a typical undulatory timescale, cycle ; for a nematode this is cycle ≈ 1 s. ¶ We can extend this idea (and thus the inertialess locomotion assumption) to granular undulating systems. We justify our extension by estimating the ratio of inertial to frictional forces in Coulomb friction-dominated systems (which are approximately rate independent):

, where the numerator is the characteristic dynamic inertial forces and m, v 0 , and cycle are body mass, average speed, and temporal period respectively; the denominator is the characteristic frictional force where and g are the friction coefficient and gravitational acceleration constants respectively. We can rewrite the above as

where the numerator can be interpreted as the time required to go from steady-state locomotion to a complete stop. Because force in a frictional fluid environment is approximately rate-independent, we have

In doing so, this ratio is then (in frictiondominated systems) exactly C. Thus, like in the viscous swimmer, in the macroscopic granular swimmers we have analyzed, C is sufficiently small (order 0.1) such that we can neglect inertial effects in granular locomotion (SI Appendix, Tables S1 and S2).

For forces acting on the body of a viscous drag-dominated nematode, we use Stokes drag from measurements from refs. (60,63). In the present study, these values likely represent an approximation of the drag forces, due to the presence of surfaceinduced hydrodynamic effects arising when swimming near a glass substrate. Regarding forces on the granular swimmers in the locomotion systems we studied, elemental forces were determined from previously developed empirical force relations (50,58), with modifications made (SI Appendix, section 4) to account for contact dynamics at granular surfaces (Fig. 3A). Knowing the elemental forces at each point along the body (depicted in Fig. 3B), the instantaneous swimming speed, ġb (t), satisfying the force balance in Eq. 2, can be numerically determined (SI Appendix, sections 4 and 5). Notably, ġb = [ x , y , ], where each element denotes an instantaneous velocity component in the forward, lateral, and rotational directions expressed in the time-varying local frame of the locomotor (which does not always align with the laboratory frame of reference). It is worth noting that the evaluation of instantaneous velocity is independent of the position of the locomotor with respect to the laboratory frame of reference (g = [x, y, ]). The net displacement over a cycle can be numerically obtained by solving the ordinary differential equation: Our previous work (49,50,55) demonstrated that RFT agrees with experimentally measured forward speeds (displayed here as body lengths per undulation cycle) of the sandfish lizard and the shovel-nosed snake. RFT calculations predict that intermediate body curvatures result in the largest displacements, and previous studies demonstrated that animals move using body curvatures that nearly coincide with the RFT-predicted speed-maximizing shapes. Given the success of RFT, we will model environmental forces with this approach.

RFT does not immediately facilitate ready understanding of how variation in paths in the configuration space results in different amounts of displacement (e.g., what happens to displacement when we go from circle paths to ellipse paths or vary parameters like radius of circle), nor does it permit rapid understanding of why certain paths could be better than others (39,45). That is, beyond merely rationalizing observed gaits, what are the principles by which we can predict paths for other living (and ultimately nonliving systems like robots) to achieve optimal performance without having to do laborious calculations? While we will not fully address these questions in this paper, we now illustrate how progress can be made using the geometric phase approach.

The first step in applying the geometric framework to our systems is to movement as a series of small displacements (translations or rotations), each introduced by small body configuration changes ("self-deformations"). We seek a mapping that relates changes in body configuration space to changes in realworld space. To construct such a map (which will prove valuable in subsequent sections as we build up the machinery of the geometric approach), we next introduce a fundamental assumption in the geometric theory. That is, we assume that the small changes in displacement (small body velocities) are linearly related to small shape changes (small "shape" velocities) via the following:

where = [w 1 , w 2 ] T is the shape of the system (Fig. 4-i); ˙ is the shape velocity, the speed with which the body curvatures are changing; and A( ) is the local connection, which encodes environmental constraint forces that relate changes in body shape to the changes in position that they induce (30)(31)(32)(33)(34)36). Each row of the local connection A( ) can the be visualized as a vector field in the shape space (Fig. 4-ii).

In our previous work on geometric methods in locomotion, we showed the validity of the linearity in Eq. 4 of a granular Purcell's three-link swimmer (39), which only has two joints (internal shape Dof = 2). Here, we consider a two-dimensional shape space identified by PCA as discussed earlier. We numerically obtained the local connection matrix A using the same approach as discussed in ref. 39. The validity of the linearity assumption for the two-dimensional reduced shape space representation is shown in SI Appendix, Figs. S1-S4 (64,65).

Given the local connection, we aim to identify relationships between the geometries of gaits and the displacements (net body translations or rotations) they produce (the accumulated geometric phase). Notably, Eq. 3 is nonlinear and time-dependent, causing additional challenge to directly calculate the geometric phase. To simplify our analysis, we can approximate Eq. 3 as

where is now time-invariant. Thus, the net displacement Δ is uniquely determined by the trajectory of shape change. The establishment of such approximation us to evaluate the geometric phase as the linear, time-independent line integral over the local connection vector field Fig. 4-iii. To further visualize and analyze the geometric phase, it can be convenient to turn the line integral into a surface integral via Stokes' theorem (38,66). Thus, the net displacement can be approximated by taking the integral of the curvature of A over the region of the shape space enclosed by a gait,

Taking the individual components of ∇ × A as height functions yields gait-independent signed scalar maps that provide an intuitive visual depiction of how shape changes relate to net motions (38,39,67). Surface integrals over the height functions predict displacements caused by the cyclic sequences of selfdeformation described by the boundary of the enclosed area. As a result, the height function provides a way to identify gaits which produce no displacements (i.e., enclose no net curvature) as well as gaits that yield large displacements (i.e., enclose significant net signed curvature).

The primary source of error in our approximation Eq. 5 is the fact that matrix multiplication is not commutative (38). This effect can be seen in the context of "parallel parking": a car cannot move sideways, but the interplay of forward and rotational velocities can generate an emergent lateral displacement if the forward and rotational velocities are properly sequenced. Hatton et al. (67) introduced the notion of systematically selecting the system gauge (the choice of body frame) such that the parallelparking effect is suppressed. The chosen gauge is generally close to, but not exactly, the center of mass and mean orientation of the body elements.

We used Eq. 6 to obtain the height functions shown in the contour plots of Fig. 5-i with average animal gaits overlaid in blue. Integrating the height function circles of different radii provides predictions of how displacements per gait cycle, Δ, depend on the maximum local curvature along the body of the animal, m s (Fig. 5-ii).

For the subsurface movement of the sandfish lizard S. scincus in a frictional fluid (Movie S1 and Fig. 2 B-ii) and for swimming of C. elegans in a true fluid (Movie S2 and Fig. 2 A-ii), comparisons of animal performance, direct RFT simulations, and height function surface integrals reveal that displacements per cycle are close to predictions (Fig. 5 A and B). In the case of the sandfish, biological data from previous studies (49,50,68) provide a range of local curvatures (SI Appendix, section 3) that are predicted to yield near-maximal displacements per cycle. This is in accord with previous muscle activation measurements that identified these template parameters as targets for the neuromechanical controller (68). In the case of the nematode, animal performance is in agreement with RFT and height function integral predictions. It is important to note here that because of power limitations in living (or synthetic) systems, displacement per cycle is not necessarily equivalent to speed. Since power generation capabilities of a swimmer are finite (e.g., muscles are not infinitely strong), larger shape changes (and therefore larger amplitude cycles) require more time to execute (69). As a result, in the case of the sandfish, the peak power-limited speed occurs at a slightly smaller amplitude. In contrast, in low viscosity regimes C. elegans is not power limited (59). Power limits arise in higher viscosity regimes (e.g., agar), where the nematode uses greater muscle power to deform both its body and the surrounding fluid (59). The shovel-nosed snake Ch. occipitalis used different waves for locomotion on the surface of sand and within sand (Fig. 2 C and D). For subsurface movement of Ch. occipitalis, the first two PCs only capture 57.6% of the variation in body configurations (Movie S3). One possible explanation for the apparent increase in postural complexity in subsurface movement may result from the increase in spatial frequency relative to surface crawling. As the number of waves along the body increases, traveling waves may become localized to such a degree that the coherence across the entire body is lost. Another possible mechanism may arise from the finite muscle torque provided by the organism in the increased stiffness subsurface environment. Torque limits may cause localized failures to realize a target waveform which may act as a random perturbation causing wave decoherence.

Despite the lack of clean circularity in the configuration space, previous work showed that RFT and experiment are in good agreement (50) which we attribute to the locality of granular resistive forces and the low slip locomotion. To test the efficacy of the geometric scheme in such a situation, we constructed the height function using the shape space spanned by the two dominant PCs, and the fully immersed environmental stresses (Fig. 3A) is shown in Fig. 5 C-i. Predicted displacements from direct RFT simulations and height function surface integrals are in agreement with animal performance (Fig. 5 C-ii and SI Appendix, sections 4 and 5) and reveal that animals use postural dynamics predicted to yield near-maximal displacement per cycle.

The surface waveform used by Ch. occipitalis, which has fewer waves and lower curvatures, produces low-slip movement that leaves behind a well-defined track of depth ≈5 mm (55) (Movie S4). Given that RFT measurements for movement at the surface used a flat plate intruder, we added an additional term to the measured RFT relations to account for the kinetic Coulomb friction drag that opposes the motion of the local segment (SI Appendix, section 4). Predictions from direct RFT simulation and height function surface integrals are in agreement with animal performance (Fig. 5 D-ii). We can rationalize the difference in the subsurface and surface waveforms by considering that granular drag force increases with intruder depth. For a fixed depth, increasing wave curvature and/or wavenumber decreases the amount of torque joints must produce (55). Thus, when moving subsurface the snake can contend with the increased environmental forces by adjusting its waveform to reduce the torque required.

To navigate in complex terrain, effective turning (reorientation combined with or without translation) is as important as translation (21,52,70). However, unlike in forward motion, both the wavelength and the amplitude of the worm body wave become time dependent (52), leading to additional challenges to reconstruct the turning behavior. It was hypothesized that additional PCs (eigenworms) are needed to fully describe the turning behaviors (52). Here, we hypothesize that turning and forward behaviors may share some common dynamics. Specifically, we posited that turning is a modulation of the serpenoid template in the shape space that can be rationalized by the geometric phase approach. Note that to fully describe worm turning behaviors on agar, we require 4 PC modes, as discussed in refs. 71 and 72. However, in this study, we demonstrate that 2 PC modes are sufficient to characterize worm turning in buffer fluid.

The nematode exhibits extraordinary maneuverability in part because of its ability to perform omega turns (73). The turning motion is called an omega turn because during the course of turning the anterior end of the body (head) sweeps near the posterior end of the body (tail), inscribing an "Omega" (Ω) shape (Fig. 6A). In Fig. 6A, we illustrate an example of an omega turn in C. elegans. It is worth noting that the worm body in Ω shape cannot be readily prescribed by a sinusoidal function. During an omega turn, the worm body orientation experiences significant (typically over 60 • ) rotation with negligible translation (typically less than 0.1 BL/cycle; see Fig. 6A). We refer to this type of nontranslating turning behavior as an in-place turn. The efficacy of such turns has made them targets for control of turning in robots (71,72). We perform PCA to explore the omega turn kinematics and notice that the modes of C. elegans forward motion and turning motion are surprisingly similar. Specifically, the two standard sinusoidal modes identified in C. elegans forward motion (Fig. 2 A-iii) can explain over 80% of the total variation of curvatures in turning motion. In Fig. 6B, we compare the PCs calculated from isolated recordings of turning motion (solid thick curves) and the standard sinusoidal basis derived from isolated recordings of forward motion (dashed thin curves, same as Fig. 2 A-iii).

In omega turn kinematics, there is a significant difference between the variance explained by the first PC (which we will refer to as mode 1) and the second PC (which we will refer to as mode 2), hence their order is not arbitrary, in contrast to the case of forward crawling described above. The similarity of the PC modes in forward and turning motions suggests a conserved shape space can realize a variety of behaviors.

Unlike in forward motion, for turning behaviors we notice that the trajectories in the shape space are centered off-origin along the w 1 axis (hence the difference in the amount of variance explained by each PC). Notably, clockwise (CW) turns are typically associated with positive w 1 offset and counterclockwise (CCW) with a negative offset (Fig. 6C). For each trial (20 individuals, one trial per individual), we calculate the distance of the center of the recorded space shape trajectory from the origin (x c ) and measure the net rotation the space. Fig. 6E shows clear correlation between x c and the body rotation (each trial measurement is represented as a blue dot).

In the 2D plane, there are three connection vector fields (corresponding to forward, lateral, and turning dynamics, Eq. 3). We can use the same methods narrated above to generate height functions for each. Performing these calculations (SI Appendix, section 5E), the rotational height function reveals two distinct positive and negative regions along the w 1 axis (Fig. 6D). To quantify the turning modulation, we perform the surface integral over standard circular templates (with a radius of m s = 8) subject to different offsets in w 1 axis. We bound the uncertainty in effective drag anisotropy by computing a series of height function with different drag anisotropy (d a ). We compare the geometric phase approach predictions on displacement and rotation with the empirically measured data in Fig. 6E and observe good agreement. Thus, a shape space made with a common set of modes can describe both worm forward swimming and turning in fluids, simply by applying the appropriate height function; near-optimal turns can be achieved with the modes from forward crawling through amplitude modulation that serves to offset circular templates gaits from the origin. This suggests that to turn, organisms may capitalize on common neural means of coordination associated with forward locomotion, as turns and forward crawling share a similar low-dimensional representation.

Thus far in the paper, systems studied are assumed to maintain continuous full-body contact with the environment during selfpropulsion. However, many animals lift limbs or body portions as they move, changing their contact state throughout a gait cycle (74,75). We therefore sought to build on our previous success in applying the geometric framework to such situations in robophysical models (76)(77)(78). Following the other examples of undulatory behavior in this work, we chose an organism that modulates environmental contact within a flowable resistive environment with vertical waves, the rattlesnake, Cr. cerastes (Fig. 7A). This organism encounters sandy substrates in its native North American desert habitat and moves by sidewinding. Sidewinders locomote on homogeneous substrates (79,80) by propagating a wave of planar body undulation coupled to an offset wave of body lifting (Fig. 7B), resulting in each body segment being cyclically lifted clear of the substrate, moved forward, placed into a nearly static contact, then lifted again, with a slight phase offset between successive segments (75, 81-83) (Movie S5). Thus, the snake generates multiple head-to-tail propagating regions of lifted movement and nearly static ground contact and moves at a nonzero angle relative to the overall headto-tail body axis (Fig. 7C) (75,77,(81)(82)(83).

Despite the apparent complexity of these movements, our previous work indicated that the self-deformation pattern of Cr. cerastes could be characterized as a template consisting of a superposition of a planar and vertical traveling wave (21), with a phase shift of ± /2 between them (Fig. 7B). The modulation (e.g., changes in the maximum amplitude) of these waves can lead to diverse behaviors (21,77,84). However, sidewinders tend to use relatively consistent horizontal waves during forward motion and are typically thought to regulate forward speed using temporal frequency changes of the wave (75,79,84).

Indeed, when we applied PCA to horizontal wave dynamics of previously collected Cr. cerastes data, we found that, across trials dynamics of the horizontal wave consists of a circular path in a two sinuous mode configuration space (Fig. 8 A and B) of a characteristic radius (and therefore maximal body curvature). We thus posit that this circle forms a control template enabling these animals to move rapidly over loose granular surfaces.

Although the vertical body dynamics have not been carefully experimentally resolved (84), they are assumed to be a traveling wave (and thus described by two modes) that sets the periodic contact pattern (Fig. 8 A and B). On level granular media (21,84), parameters describing the vertical template remained approximately constant. Therefore, to model the vertical wave interaction, as in refs. 21, 77, and 85, we introduced a weighting prefactor, c, that specified how much of the environmental force each infinitesimal segment experienced. # Specifically, we modify the resistive force balance in Eq. 2 as

Previous work (21) revealed that the three-dimensional pose of Cr. cerastes could be represented by a horizontal wave (characterized by w 1 , w 2 ) coupled to a phase-shifted vertical wave that sets the environmental contact condition. To properly couple the contact function to the in-plane shape, we introduced the vertical wave:

where a is the amplitude of the vertical wave, w 1 and w 2 describe the in-plane wave shape. To set the contact using the vertical wave description, , we defined the smoothly varying function c(s) = 1/ 1 + exp[ (s) + b] , where c ∈ [0, 1] sets the local fraction of the environmental force experienced as a function of position along the body, and b sets contact width. To be consistent with previous observations, a = 15 and b = 0.5 are chosen so that approximately 34% of the animal's body is on the ground (84). Fig. 7B shows how environmental contact couples to an in-plane shape, and Fig. 8C shows how this contact varies throughout the in-plane shape space for an animal with n = 1.5 waves along its body. Fig. 8D shows four RFT simulation snapshots throughout one undulation cycle. Contact patches (dark regions) originate near the head and are propagated toward the tail. Given the experimentally observed oblique direction of travel (relative to the head-to-tail body axis), we expect the kinematics in our modeling to produce significant displacements in both the x (forward) and y (lateral) directions. We therefore numerically computed connection vector fields (Fig. 9A) and height functions to visually and intuitively prescribe motions along both the x and y directions (Fig. 9B). We define the total predicted displacement is given as Δ = (Δ 2

x +Δ 2 y ) 1/2 , where Δ x and Δ y are displacements predicted from x and y height functions, respectively. Fig. 9C shows that, for movement on granular media, direct RFT simulations (dashed tan curve) and geometric computations (solid tan curve) predict similar maximal displacements. The RFT gait predicted to yield peak performance differ slightly from those predicted to from the "cartoon model" presumably because at large amplitudes slip in the primary direction of motion occurs. Despite the differences in predicted gait amplitude, the displacement curves predicted are not highly sensitive to gait amplitude variation over a broad range. Note that the discrepancy between RFT simulation and geometric phase approach for sidewinding on sand can be a result of the noncommutativity of body velocities. As shown in Fig. 9A, the body velocity in xand y-directions have comparable magnitudes, which can lead to relatively large noncommutativity effect in body velocities (38).

Differential Turns in Sidewinding. As with the worms, the geometric phase approach can also help rationalize the spectrum of "differential turns" observed in Cr. cerastes (21). Such turns are interesting because those sidewinders can modulate the net translational displacement associated with a particular turning angle. That is, sharp differential turns (e.g., ∼90 • per cycle) are often accompanied by reduced translational displacement, and gradual differential turns (e.g., less than 20 • per cycle) by large translation. Therefore, unlike in straight sidewinding where the animals use relatively consistent horizontal waves (e.g., consistent wave amplitude and propagation speed), animals exhibit varying horizontal wave dynamics during differential turns.

The differential turning mode, first analyzed by Astley et al. horizontal wave (Fig. 10A). Depending on the of modulation, can control the degree of turning during translation. Notably, in prior work, the amplitude modulation refers to the modulation of elemental velocity distribution (from anterior to posterior) (21). It is yet not clear how internal curvature should be adapted to facilitate such amplitude modulation on spatial velocity distribution.

We used PCA to find the curvature modes during differential turns. We noticed that the first two principal components can account for over 69.1% of the variance (Fig. 10B); and the two modes for straight sidewinding are almost identical to those for differential turns with two subtle differences: i) the spatial frequency changes from 1.5 in straight sidewinding 1.2 in differential turn, and ii) the phasing relationship between two modes changes from mode 1 ahead of mode 2 in straight sidewinding to mode 2 ahead of mode 1 in differential turn. We suspect that these subtle changes emerged from the variation in the horizontal wave. We project body curvatures during differential turns onto the first two PC modes (Fig. 10C). We notice that the configuration space trajectory is approximately circular, with its center offset from the origin. As with the worm turning, we posit that the offset of the center from the origin can serve as an indicator of the degree of turning and translation in the snakes. We measure the rotation and translation for each cycle (95 cycles over 47 trials, a trial might include multiple cycles). The translation (Fig. 10 E, Left) and rotation (Fig. 10 E, Right) in each cycle (represented as a blue dot) are plotted as a function of w 1 -offset of the trajectory center (arithmetic average) from the origin. Linear regression between w 1 -offset and rotation shows significant relationships (translation: r 2 = 0.20, P < 0.01; rotation: r 2 = 0.50, P < 0.01).

We use the geometric phase approach to explain the observed correlation. We numerically compute the height functions (Fig. 10D) on granular media using the same contact function and RFT relationships as straight sidewinding. The clusters of positive and negative volumes are distributed along the axis of w 2 = 0 in height function, indicating that the introduction of offset in w 1 direction can indeed lead to body rotation.

We then perform surface integrals over circular gaits on the height functions. Specifically, we used the following equations to prescribe off-centered circles in the calculation:

[8]

We compare the surface integral with RFT calculation (integrating Eq. 8 over a cycle t ∈ [0 2 )) and the fitted linear regression, and observed good agreement. || Our analysis illustrates that the seemingly distinct behaviors of straight sidewinding and differential turn (through amplitude modulation on spatial velocity distribution) share the same shape space (PC modes). Moreover, the simple modulation scheme (with a single variable: w 1 -offset) reconstructs the complicated spectrum of behaviors with different degrees of displacement and rotation, which further facilitates a relatively simple understanding of seemingly complex sidewinder locomotion.

In this paper, we presented a biological experimental application of a geometric framework of locomotion proposed in the 1980s by physicists (25) to describe movement and developed during the last three decades for robotic applications by control theorists (29,31,38). Application of this framework to organisms across scales and levels of complexity revealed that selfdeformation kinematics were well described by serpenoid templates (near circular paths in low-dimensional body configuration spaces). Further, observed animal self-deformation kinematics and locomotor performance coincided with predictions that nearly maximized the surface integral over this curve in a diagram called a height function, which corresponds to nearly maximizing a geometric phase in the space of animal body configurations. We can thus posit that a emergent guiding principle for control of limbless undulatory locomotion in highly damped environments || Note that forward and lateral height functions in Fig. 10D have opposite sign to those in Fig. 9 because of the changes in the relative phasing relationship between mode 1 and 2.

is: maximize geometric phase relative to the time/effort needed to complete the cycle, using a template which is approximately a circle in a two-dimensional configuration space. This is true both in situations with continuous and variable environmental contact-for cm-scale lizards and snakes in flowable (granular) and frictional environments (both of which are ubiquitous in natural organism environments) as well as a tiny nematode worm in a fluid environment. Finally, modulations to serpenoid templates can explain the turning behaviors (omega turn and differential turn), which further indicates the generality our templates of the underlying physiological, neural, and biomechanical mechanisms responsible for the selfdeformation patterns.

Why is this seemingly abstract approach, requiring mathematical tools not traditionally represented in the fields of bio-and neuromechanics, of value in locomotion analysis? First, recognition that the dimensionality reduction scheme proposed in ref. 52 applies to diverse and significantly more complex organisms than nematode worms is a useful step in discovery of candidate templates (i.e., high-level control targets). Second, the diagrammatic approach simplifies search for candidate paths in configuration spaces; typically such a search requires brute force computation of all allowed paths but with height functions in hand, optimal selfdeformation patterns for translation rotation (71)] becomes relatively straightforward to hypothesize. Third, the scheme can be readily adapted (69) to imposing biologically relevant constraints like internal force, power, or energetic limitations.

Amplifying on this last point, to demonstrate our method, in this paper, all computations were performed in the space of shape changes, where all shape changes are equally easy to execute and the geometric phase nearly maximized displacement per cycle. But recent work (69) has demonstrated connection vector fields can be computed within modified metric spaces (such as weighting shape changes by internal force and power requirements). Modifying the underlying metric to account for physiologically relevant constraints could provide insights into the goals and limitations associated with a broader range of behavior and inherent biological variability. With such theoretical tools in concert with new experimental tools, we can now address how lower-level physiological, neural, and biomechanical mechanisms ["anchors" in the parlance of Full and Koditschek (13)] conspire to generate the template dynamics (18,68). In large-scale organisms, electromyographic and neural recording tools at the macroscale have been used for decades to assay neuromechanical control across taxa (86). Relevant to animals studied here, such tools proved useful to test a hypothesis of shape control in sandfish locomotion (68). However, such tools provide crude assays relative to those available in model organisms (like C. elegans). The optical transparency and genetic mutability of these worms provide experimental opportunities to connect templates to underlying anchors through genetic circuit manipulation, optogenetics, and calcium imaging. Using the geometric scheme to develop templates for diverse environments, we can, for example, test the hypothesis that nematodes control for force and thus shapes emergently arise from neuromechanical feedback (10,87).

Finally, the results in this paper used reductions of the kinematics of animal movements to simple patterns of self-deformation and used RFT to describe environmental interactions. We posit that this is because we focused on relatively simple tasks, such as escape and steady transit in homogeneous media. The extent to which more complex locomotor behaviors (e.g., refs. 88 and 89) could be amenable to such simplification remains an open question. Complexity may arise from unusual gait paths in the shape spaces or from additional modes. While C. elegans uses a circular gait in agar (52), changes in the viscosity of the environment could lead to a spectrum of gaits, ranging from circular to elliptical gaits. Such gait spectrum is observed in the continuum of lizard body elongation and limb reduction, where elliptical gaits emerge as limb size reduces (22). Generalizations to more complex descriptions (e.g., more modes) are straightforwardgeometric methods can handle higher dimension although the visualization of height functions becomes more difficult (90,91). It will be interesting to explore other, potentially more complex, locomotor behaviors not described by planar traveling waves of body bends [e.g., rectilinear motion in snakes and peristalsis in worms, walking and crutching of mudskippers (76)].

For all experiments, high-speed videos (excluding C. elegans) were recorded as animals moved on various substrates, and animal postures were obtained by identifying tracking features located along the midline of the dorsal side (lateral side for C. elegans) of the animal in each frame. For subsurface animals, lead markers glued to the animal were visible in X-ray images and extracted as a part of a previous study (50). For surface animals, either markers or features on the snakes were tracked through time, and again these points were interpolated using a cubic spline. For C. elegans, animals were placed in 10 μL of S-basal on a glass slide. 3D movement was constrained by a glass coverslip to a 50 μm height, with the use of tape (Kapton). Worms were imaged using a brightfield microscope (Leica ATC 2000). For forward crawling data, midlines were obtained using custom MATLAB code. Binary masks of the organism were created via thresholding. The binary masks were then skeletonized and splined. For turning data, DeepLabCut (92) was used to track the self-occluding postures where binarization and skeletonization fails. For Ch. occipitalis on the surface of sand, the positions of naturally occurring evenly spaced black bands were identified in each video image frame as part of a previous study (55). For Cr. cerastes, infraredreflective markers were placed along the animal, and a Natural Point Optitrack Flex 13 camera system automatically identified and recorded marker positions at 120 frames per second. Further biological information and experimental details for S. scincus and Ch. occipitalis are given in SI Appendix, Table S1 and for Cr. cerastes in SI Appendix, Table S2. Details on numerical analysis and geometric phase approach are provided in SI Appendix, sections 3-5.