A physical theory of larval Drosophila behaviour

All animal movement must ultimately be governed by physical laws. As a basis for understanding the interactions between the nervous system, musculature, body mechanics and the environment that govern behaviour in the fruit fly larva, we here develop an effective theory for the physics of its motion in three dimensions. We start by defining a set of fields which quantify stretching, bending and twisting along the larva’s anteroposterior (AP) axis. We then perform a search in the space of possible physical theories that could govern these fields, by using symmetry considerations, stability requirements and physical reasoning to rule out possible terms in our theory’s Lagrangian (which governs its energy-conservative physics) and Rayleigh function (which governs its energy-dissipative physics). We restrict attention to the physics that dominates at long-wavelengths, which allows us to arrive at a unique, simple theory of the larval midline, governed by a minimum of phenomenological parameters that capture both purely biomechanical as well as neuromuscular effects. Owing to the simplicity of our theory, we are able to derive most of our results analytically. The model makes strong quantitative predictions for the dynamics of peristalsis, rolling, and self-righting, and also successfully predicts statistical properties of these behaviours and of unbiased substrate exploration.


Introduction
Animals have evolved an incredibly diverse set of strategies for locomotion. This allows them to search for resources, adequate habitats, or mating partners, as well as escape predators. Yet, despite this diversity, all animal movement is fundamentally a physical act, and must therefore be governed by physical laws.
Determining the exact form of these laws, which must arise through the interactions between the nervous system, musculature, body mechanics, and the environment [1,2], is a fundamental and inspiring problem in modern biology.
Here we develop an effective theory [3] describing the physics of 3-dimensional behaviour in larvae of the fruit fly, Drosophila melanogaster. An effective theory aims to provide an appropriate, abstracted description of the important physics governing a system's behaviour, in the form of a computationally tractable model that advances intuitive understanding at a relevant scale, while remaining neutral on the causal counterparts of the mechanisms governing the results. This approach, widely adopted in physics, has the advantage of facilitating both concrete predictions that can be compared to measurements on a specific system (with minimal parameter fitting) and substantial generality of application and insight.
One aim of this paper is to show that it is possible to develop such a theory for the complex dynamics of behaviour of a whole organism. We take as our target system the Drosophila larva. This animal hatches and grows within decaying fruit before travelling further afield in order to pupate and metamorphosize into its adult form, and so must be capable of moving in complex three dimensional environments. Correspondingly, it has been observed to perform a wide repertoire of behaviours under experimental conditions [4], including peristaltic crawling [5], rolling [6][7][8], self-righting [9][10][11], rearing [4], hunching [7], and digging [12][13][14][15].
As well as extensive behavioural data, there is currently a concerted effort to map the neural connectome of the larva, and to link observed neural anatomy and dynamics to behaviour [16][17][18][19][20][21][22][23][24][25][26][27]. Moreover, recent modelling studies have already explored the importance of biomechanics in understanding the animal's crawling (in 1 spatial dimension) [28][29][30][31] and substrate exploration/taxis behaviours (in 2 spatial dimensions) [30,32,33], emphasising the crucial role that sensory feedback, which conveys some measurement of the physical state of the animal to its nervous system, plays during behaviour [28][29][30][31]33]. Cutting-edge research is beginning to experimentally probe the physics of the animal's body [31]. The larva thus offers a number of advantages as a model system for studying the physics of behaviour.
Additionally, the larva provides a prototypical example of the bilaterally symmetric, slender, soft, segmented body physics that are common to very many animals. The larva's body is entirely soft, supported by a "hydrostatic skeleton" consisting of a continuous fluid-filled cavity that runs through the body, and which mechanically couples the head and tail during peristaltic crawling (acting as a "visceral piston" [5], similar to that observed in caterpillars [34,35]). The internal fluid is enclosed by segmentally repeating outer tissue layers containing the musculature and cuticle. A theory of the soft body physics of the larva should thus advance general understanding of soft physics at organismal scales, and also has engineering potential, for instance in the emerging field of soft robotics.
In outline, our aim is to construct a theory to describe the motion of the larva's anteroposterior (AP) axis (which will we refer to as the larval 'midline') in 3-dimensional space, including torsional motions around this axis. Our approach hinges on the principled construction of a dissipative Lagrangian field theory that has a minimum of phenomenological parameters, and is able to capture both purely biomechanical as well as neuromuscular effects. We derive quadratic ("free-field", or non-interacting) Lagrangian and Rayleigh dissipation functions, which describe the storage of energy or the flow of power, respectively, and which give rise to simple linear dynamics. Owing to the simplicity of our theory, we are able to derive almost all of our results analytically, without having to resort to simulation or numerics. We show that our effective theory is able to explain and predict many features across a broad range of larval behaviours, including selfrighting, hunching, rearing, peristalsis, rolling, and unbiased substrate exploration. We find evidence that the statistics of several real larval behaviours can be surprisingly well described by equilibrium statistical mechanics, with elementary low-energy, long-wavelength motions seeming to dominate the animal's behaviour due to Boltzmann suppression of the higher-energy, shorter-wavelength motions. We discuss the implications of these results for gaining insight into behavioural variety and control, both in larva and in other animals.
2 Experimental methodology 2.1 Self-righting behaviour assays Self-righting (SR) assays were conducted as done previously [9] on 1st instar wild-type larvae (w1118), separately to our experiments on unbiased behaviour and rolling. Briefly, parental lines were raised at 25 • C in collection cages bearing apple juice-based medium agar plates, supplemented with yeast paste. From these plates, stage 17 embryos [36] were collected and transferred to a fresh plate. Next, freshly-hatched first instar larvae were transferred again on 0.9% agarose plates. After one minute of acclimatisation, we conducted self-righting assays by gently rolling over the larvae with a small brush pen to an upside-down position. Wild type larvae normally take approximately 6-7 seconds to rectify their position (self-right). SR sequences were recorded at 10 fps with a Basler ace acA800-510µm (CMOS) monochrome USB 3.0 Camera mounted on the Leica MZ75 stereomicroscope.

Peristalsis assay
The data we used to quantify peristalsis behaviour was previously collected and published in [37]. Videos of forward and backward crawling were collected at 30 fps. The lengths of abdominal segments were measured at the left and right sides of the body; we used the mean of left and right measurements to estimate the midline segment length prior to further analysis (see below).

Fly strains
For behavioural assay of unbiased behaviour and rolling, we used R69F06-GAL4 and a control line with no GAL4 expression, w;;attp2, from the Rubin GAL4/LexA collection [38,39]. These were crossed to UAS-CsChrimson::mVenus [40] for optogenetic activation of rolling behaviour. For mechanical nociceptive stimulation behavior experiment, CantonS were used. All flies are raised at 25 • C for 4 days -3rd instar larvae before wandering stage animals are used for behaviour experiment.

Behavioural apparatus
The apparatus comprises a video camera (FLIR GS3-U3-51S5M-C camera, 2048 × 2048), for monitoring larvae, a light illuminator (LED 624 nm) for optogenetic activation, and a computer, similar to [7]. Recordings were captured at 30 fps and were controlled through the Multi-Worm Tracker (MWT) software (http: //sourceforge.net/projects/mwt) [41], whilst control of the hardware module was controlled through the Stimulus Control Module (SCM) software. For mechanical stimulation, video was captured 30 fps using FLIR software.

Behavioural experiments
Embryos were collected for 24 hours at 25 • C. Foraging third instar larvae were used for all experiments. Larvae were raised in the dark at 25 • C for 4 days on fly food containing trans-retinal (Sigma, R2500), at a concentration of 200 µM. Before an experiment, larvae were separated from food by suspension in 15% sucrose and then washed with water. Larvae were dried, then transferred to the centre of a 25 × 25 cm transparent plastic, square plate covered in a layer of 2% agar gel. Between 20-60 larvae were transferred to the plate for any given recording. Optogenetic stimulation was delivered in two 15 second bouts at an irradiance of 600 µW/cm2, with a 30 second interval between bouts. Experiments were run under infrared light to avoid tonic activation of UAS-CsChrimson::mVenus. For mechanical nociceptive stimulation, single larvae were placed on 2% agar gel. After larvae started crawling, they were pinched by forceps (Dumont, N5) to evoke rolling behavior.

Behaviour quantification
For calculation of larval midlines, behavioural recordings were captured with the Multi-worm Tracker (MWT) software http://sourceforge.net/projects/mwt [41]. MWT returns a 2D contour (outline) for each larva tracked, at 30fps. From these contours, we computed the larval midline using the choreography package (bundled with MWT) [41]. The estimated midline is given by a set of equally spaced points, allowing characterisation of transverse bending but not axial compression/expansion. Larval midlines and contours were used to score bouts of rolling behaviour, as described elsewhere [27]. Larvae that were tracked for fewer than 5 seconds, or travelled less than one body length in distance, were rejected. To determine angular velocity of rolling behaviour, videos were manually scored using the Fiji software (NIH) [42]. First, we scored the angle of larval rotation to the nearest 60 degree interval (0, 60, 120, 180, 240, 300), where 0 degrees represented dorsal-side up, and 180 degrees represented dorsal-side down. The trachea and posterior spiracles were used as landmarks to determine the degree of rotation around the anteroposterior axis. The duration of a roll was demarcated by the first and last frame in which the larva presented movement perpendicular to it's body axis. The mean angular velocity of rolling was calculated as the total angle rotated divided by the duration of roll.

Data analysis
For unbiased behaviour and rolling behaviour, postural principal components ("eigenmaggots") were extracted from the midline data using the singular value decomposition (SVD)-based principal components analysis (PCA) provided by the scikit-learn python machine learning module [43]. Prior to PCA, we removed the overall translational and rotational degrees of freedom of the animal by computing the angles between segments [44].
For peristalsis behaviour, principal components were extracted from measurements of abdominal segment lengths during forward and backward crawling (data previously published in [37]) using the same SVD-based PCA decomposition as for unbiased behaviour and rolling behaviour. Overall translation was not measured in this dataset so did not need to be removed prior to PCA.
For comparison to our theory, the experimental data was also projected onto the modal basis of our effective theory. We projected the abdominal segment length data onto the eigenvectors of the circulant second difference matrix D 2,c , which is identical to a spatial discrete Fourier transform basis (Appendix E), to obtain estimates of the axial modal coordinates X i . We normalised to correct for the truncation of the axial modal basis (since only abdominal segments were measured), and for figure 6 applied a brickwall (0.5-1.5Hz) bandpass filter to the first pair of axial modes to remove noise and artifacts associated with the truncation. Since we have measurements for seven abdominal segments we are able to project only onto the first six non-zero-frequency modes. Higher modes than this become indistinguishable from lower modes due to the reduced number of samples available and the Nyquist sampling theorem.
The midline data from unbiased behaviour and rolling behaviour experiments was projected onto the eigenvectors of the free-boundary fourth difference matrix D 4,f . We estimated the eigenvectors numerically using the numpy python numerics module [45] as the resulting eigenvectors are more accurate than our analytical approximation. To remove the overall translational and rotational degrees of freedom we performed this projection using the angles between midline segments, as for PCA (above), by similarly transforming the eigenvectors into the angle space. Unbiased behaviour occurs mainly in the plane of the substrate with little rolling, so we use the projection to estimate the mediolateral transverse modal coordinates Y i during this behaviour. By contrast, the larva continuously rotates during rolling so that our top-down view allows us to effectively estimate the modal amplitude in the mediolateral-dorsoventral plane, To estimate the mean modal frequencies during behaviour we first estimated the power spectral density (PSD) of each modal coordinate, for each individual larva tested, using the Welch spectrogram averaging method, implemented in the scipy python scientific computing module [46]. Since peristalsis is a highly rhythmic behaviour, the estimated PSDs are tightly peaked [37]. We therefore estimated the mean modal frequencies during this behaviour using the strongest peak in the PSD. Unbiased behaviour is aperiodic, displaying a "spread" in the PSD estimates rather than cleanly localised peaks [32]. We therefore estimated the mean modal frequencies during this behaviour using the spectral centroid, i.e. the weighted average of the PSD. The higher modes are associated with a very low variance, and for the fifth transverse mode and above the spectrum drops towards the uniform PSD expected of white noise. In this case the spectral centroid stays close to the mid-point of the spectrum; for our sampling rate of 30 Hz, the Nyquist frequency is 15 Hz and the spectral centroid of white noise is 7.5 Hz. This explains the plateau in the measured mean frequency of the higher modes at around this frequency.
Unless otherwise noted above all statistical analyses and tests were performed using the scipy statistics submodule [46].

An effective theory of the larval midline
In this paper we present a simple, effective physical theory of larval behaviour that captures the effects of the passive mechanics of the larva's body as well as it's neuromuscular system. In this section we will provide an accessible summary of the construction of our theory from first principles, while a detailed derivation is provided in Appendix A. We provide some intuition as to why it is possible for neuromuscular effects to be incorporated into our purely physical theory in the Discussion, and explicitly show how neuromuscular effects can be accounted for in Appendix B.
Our theory is constructed from first principles, by borrowing an approach common to several areas of modern physics including statistical mechanics [49], quantum field theory [50,51], and condensed matter [52]. Essentially, this involves a "search" in the space of possible physical theories of the midline. The key idea is to write down the most general possible theory, and then eliminate parts of it according to basic physical principles or rudimentary observations or assumptions about the system at hand. This is accomplished by first writing very general Lagrangian and Rayleigh densities. The Lagrangian density governs the conservative physics, which involves lossless storage and transformation of energy, while the Rayleigh density governs the dissipative physics, which involves the flow of power into or out of the system. The general Lagrangian and Rayleigh densities are then written in the form of an infinite sum of different terms, and we proceed to eliminate terms according to basic physical principles or rudimentary observations or assumptions about the system at hand, such as the presence of symmetries or the importance of long-wavelength behaviour. This process provides specific Lagrangian and Rayleigh densities for the larval midline.
The great benefits of this effective theory approach are that it tends to yield minimal theories governed by just a few phenomenological parameters, and it makes clear the irreducible set of assumptions governing s #(s, t) Constructing an effective theory of larval midline physics. A) An illustration of the Drosophila larva with its midline overlaid (black, blue=head, orange=tail). B) We measure the displacement of the larva from an undeformed configuration (bottom) to a deformed configuration (top) via four scalar fields parametrised by arc-length s (measured along the body from s = 0 at the tail, to s = 1 at the head) and time t. The fields correspond to the 3 translations x(s, t), y(s, t), z(s, t) and the torsional rotation of the body about each point θ(s, t). C) We restrict attention to low energy, long wavelength, low frequency motions of the midline (top) rather than the higher energy, short wavelength, high frequency motions (bottom). D) We impose continuous symmetry of the effective theory under overall translations and rotations. E) We impose discrete symmetry of the effective theory under anterior-posterior (A/P) axial reflections, left-right (L/R) mediolateral reflections, dorsal-ventral (D/V) reflections, and clockwise-counterclockwise (CW/CCW) torsional reflections. We also require that our effective theory be described by an analytic dissipative Lagrangian field theory, and ensure that the theory does not suffer from the Ostragradsky instability [47,48] (not illustrated).
the theory. Furthermore, the limitations of an effective theory tend to be fairly explicit and there is a clear way to systematically extend the theory through re-incorporation of additional terms that were initially ruled out [3]. The construction of such effective physical theories has been extremely successful in improving our understanding of biophysics at a molecular scale [53,54], for instance in predicting the physics of long biopolymers such as DNA without having to take into account the detail of electrochemical interactions, etc; here we hope to show that the approach can be successful also at a whole-organism scale.
Here, we we will summarise the steps taken in our derivation and highlight the assumptions that are needed at each step. Our approach proceeds as follows: • We introduce four scalar fields characterising the axial motion, mediolateral/dorsoventral transverse motion, and torsional motion of the midline relative to an undeformed configuration ( Figure 1A). We assume that the full deformation of the midline can be characterised in terms of these four fields and the midline is "straightened out" and "untwisted" in its undeformed configuration.
• We write the most general possible physical theory of the midline by introducing a Lagrangian density and a Rayleigh dissipation density that depend upon the fields and their derivatives with respect to space and time. The Lagrangian density characterises the energy-conservative physics of the midline while the Rayleigh density governs the energy-dissipative physics. We thus assume that the physics of the midline can be formulated as a dissipative Lagrangian field theory.
• We expand the Lagrangian and Rayleigh densities as Taylor series in the fields and their derivatives, assuming that the Lagrangian and Rayleigh densities are analytic. This step gives us a concrete representation of the theory with which to work, by providing us with an exhaustive set of terms that can be included into a general theory of the midline.
• We split our Lagrangian density into four separate Lagrangian densities governing the fields alone, along with an "interaction" Lagrangian density governing interactions between different fields. We perform the same decomposition for the Rayleigh density. This decomposition requires no additional assumptions, it is simply a convenient way to arrange the terms of the Taylor expansion. In Appendix A, Table 1 we show all 27 terms in each of the field-specific Lagrangian densities up to the fourth order.
• Focusing on the field-specific Lagrangians, we eliminate all terms including derivatives above the first order in time; these terms would lead to the Ostragradsky instability [47,48] and our theory would be non-physical.
• We then eliminate all terms violating Galilean symmetry ( Figure 1C). Retaining these terms would lead to the physics of the midline depending upon its absolute position and orientation in space. Since we are considering the physics of the midline in isolation and not considering environmental interactions, these terms must be discarded. Here our key assumption is the isotropy of space [55], which is a common assumption throughout physics and reflects our belief that there are no privileged frames of reference.
• We remove all terms violating reflection symmetries ( Figure 1D). This represents our assumption that left-and right-handed (L/R) mediolateral deflections of the midline should have equivalent physics, as should dorsal and ventral (D/V) deflections and clockwise and anticlockwise (CW/CCW) torsions. We also incorporate anterior-posterior (A/P) reflection symmetry. Note that, when combined with Galilean symmetry, these symmetries are not mutually exclusive -A/P symmetry can be seen as a consequence of mediolateral rotational and L/R reflection symmetry, or dorsoventral rotational and D/V reflection symmetry, for instance, since rotating the larva through 180 • and then flipping its "handedness" is equivalent to an A/P reflection.
• We require that the kinetic energy terms take their standard form as the square of the first time derivative of the field. This removes relativistic corrections to the kinetic energy, which are irrelevant given that the larva moves very slowly relative to the speed of light, and additionally removes a term from each of the transverse free-field Lagrangian densities, corresponding to the Rayleigh correction to the rotational inertia for beams of finite cross-section (not to be confused with the Rayleigh dissipation function) [56].
• We introduce characteristic length and time scales and retain only those terms in the expansion of lowest order in the inverse of these scales. This reflects our focus upon long-wavelength, low-frequency, low-energy motions of the body ( Figure 1B), and is motivated by the observation of considerable spatial correlation lengths in the behaviour of the larva (Figure 9, axial strain correlation length = 0.33 body lengths during forward crawling, 0.38 body lengths during backward crawling; segment orientation correlation length = 7.66 body lengths during unbiased behaviour, 0.73 body lengths during rolling).
We are left with second-order Lagrangian densities for the axial and torsional fields, and fourth-order Lagrangian densities for the transverse fields, each each including just two terms -one kinetic energy term, one potential energy term. The axial and torsional Lagrangian densities correspond to the Lagrangian of the classical wave equation while the transverse Lagrangian densities correspond to the Euler-Bernoulli beam theory.
• We look for the lowest order terms in the interaction Lagrangian density and find that all interactions must be at higher order than the long-wavelength free-field Lagrangians, i.e. interactions only occur above second order for axial and torsional fields and above fourth order for dorsoventral and mediolateral transverse fields.
• We apply the same procedure to the Rayleigh density, finding again that no interactions can occur at or below the order of the field Rayleigh densities; at long-wavelengths the non-interacting Rayleigh densities correspond to Kelvin damping.
Following these steps we are left with a long-wavelength, free-field effective theory of the larval midline in which the conservative mechanics is governed by the linear theory of elasticity [57], with axial and torsional motions obeying the classical wave equation and transverse motions obeying the Euler-Bernoulli beam theory. All four fields are subjected to linear Kelvin damping.
We note that our theory is strongly reminiscent of the phenomenological elasticity models that have been used with great success in understanding mechanical properties of DNA and other polymers [49,53,54,58], with the key difference being the inclusion of kinetic energy, owing to the larva's relatively large size and our interest in the larva's deterministic dynamics.
Since behavioural data by necessity tracks only discrete points along the larval body, we discretise the scalar fields and their derivatives using standard finite difference schemes. We choose to measure the fields at N = 12 discrete points along the midline, corresponding to the boundaries between the 11 body segments T1 -T3, A1 -A8 along with the head and tail extremities. It is necessary to introduce boundary conditions on our scalar fields in order to correctly specify the particular difference operators used in our discretisation. We impose vanishing force/moment (Neumann) boundary conditions on the transverse and torsional fields, and periodic boundary conditions on the axial field. The Neumann boundary conditions represent the absence of bending or twisting moments at the head and tail extremities, while the periodic axial boundary conditions represent the presence of a visceral piston in the larva which couples the motion of the head and tail. We shall see that several features of larval behaviour are tied to the properties of the difference operators in our effective theory.
We reduce the number of free parameters in our model by assuming that the transverse fields have an "internal" circular rotational symmetry, i.e. we assume that the effective stiffness and damping associated with mediolateral and dorsoventral motions should be equivalent. This can be motivated by the approximately circular cross-section of the larva, which should lead to an identical passive bending response in the mediolateral and dorsoventral directions. The introduction of periodic boundary conditions to model the visceral piston of the larva introduces a similar "rotational" symmetry into the effective theory. We discuss the implications of these symmetries later in the paper, during our modal analysis and while studying peristalsis and rolling behaviours.
The transverse circular symmetry removes two parameters from our theory, so that ultimately we are left with just six free parameters, three of which (the natural frequencies) govern the effective ratio of elastic to inertial energy storage, while the remaining three (the damping factors) govern the effective ratio of damping to the combined elastic and inertial energy storage. We stress that these parameters are phenomenological -for instance the natural frequencies depend upon apparent or emergent elastic properties of the body that may depend not only on the "inert" mechanical properties of the body but also upon actively controlled muscle stiffness, detailed anatomy such as the geometry of folds in the cuticle, neuromechanical entrainment, etc.
Moving on, as a final step in the construction of our theory, we convert from a Lagrangian description to a Hamiltonian description. These approaches are formally equivalent, yet have different strengths when constructing an effective theory (e.g. see discussion in [59]). On a basic, pragmatic level, the Lagrangian is specified in terms of coordinates and their time derivatives, the generalised velocities. The straightforward relationship between coordinates and generalised velocities makes it easier to determine the transformation properties of terms in the Lagrangian expansion, making it useful for the construction of our theory. By contrast, the Hamiltonian is specified in terms of coordinates and their canonically conjugate momenta. In a sense these quantities have a "life of their own" -we gain tremendous freedom in choosing variables to describe our system, simply because the momenta do not have to be the time derivatives of the coordinates. This generally makes the Hamiltonian approach useful for analysing systems, but the additional freedom can make it less useful when constructing a theory from the ground up. The Hamiltonian approach can also make it easier to leverage the conserved quantities associated with a theory, and facilitates the application of statistical mechanics.
Ultimately, we write the discretised free-field Hamiltonian of our effective theory as which corresponds to the total mechanical energy of the midline, and is the sum of kinetic and potential energies. The discretised free-field Rayleigh dissipation function of our effective theory is In these expressions the lower-case bold symbols x, y, z, θ denote vectors containing measurements of the scalar fields x, y, z, θ at discrete points in space, while the vectors p x , p y , p z , L contain the discrete measurements of the translational momentum and angular momentum. The parameters ω s , ω b , ω t measure the ratio of elastic to inertial effective forces and determine the absolute frequencies of stretching, bending, and twisting motions, respectively. Meanwhile, ζ s , ζ b , ζ t characterise the overall effective damping of these motions relative to their combined inertial and elastic energy storage. The upper-case bold symbols D 2,c , D 4,f , D 4,f , and D 2,f , denote the finite difference operators we mentioned above; we have absorbed the additional parameters associated with the discretisation into the relevant effective elasticity parameters. D 2,c is the (N − 1) × (N − 1) second difference matrix with periodic boundary conditions D 2,f is the N × N second difference matrix with free-free boundary conditions and D 4 is the N × N fourth difference matrix with free-free boundary conditions The dynamics of our theory is given by the dissipative Hamilton's equations - where q i , p i are a canonically conjugate pair of coordinates and momenta.

Modal analysis
Having formulated our effective theory, we next seek a convenient coordinate system in which to explore its behaviour. To accomplish this aim, we apply a technique called modal analysis that is similar to application of principal components analysis (PCA) to measured data. Modal analysis gives us a description of our theory in terms of a set of non-interacting, elementary collective motions of the body segments, with each collective motion measured by a modal coordinate, analogous to a principal component, and possessing a characteristic frequency and damping factor as well as a characteristic mode shape describing the spatial pattern of movement, analogous to a principal component vector. We note that this approach is quite generally applicable to mechanical systems [55], and as such is used across a wide range of disciplines and scales, for instance it is used to determine the spectral properties of small chemical compounds [44], to extract the dominant types conformational changes in proteins [60], to ascertain the vibrations of musical instruments [61], and to determine the structural resonances of buildings in civil engineering [62]. The equivalence between principal components analysis and modal analysis becomes exact for systems in thermal equilibrium, as we will explore below.
In our effective theory the mode shapes for each field and the relative relationships between their frequency and damping parameters can be determined entirely from the difference operators in the discretised freefield Hamiltonian (Equation 1). Thus, these features are completely parameter-independent. Meanwhile, the absolute frequency and damping scales for each field have a simple relationship to the parameters in the effective theory. To make this more explicit, the mode shapes are formally given by the eigenvectors of the difference operators in our effective theory, so that transforming to modal coordinates diagonalises these operators. This means that the total mechanical energy can be written as while the Rayleigh dissipation function is also diagonalised, and may be given as Here Q i denotes the i'th modal coordinate, P i denotes its conjugate momentum, ω i is the characteristic frequency of this mode, and ζ i is its damping factor. We write these parameters in terms of an absolute frequency scale ω, absolute damping scale ζ, and the i'th eigenvalue λ i of the relevant difference operator. Clearly, the ratios of characteristic frequencies and damping factors depends only on the eigenvalues of the difference operators, since we have Modal decomposition of the effective theory. A) Mode shapes predicted by our effective theory (coloured lines) closely match eigenmaggot shapes extracted from real forward/backward peristalsis (black lines, axial), unbiased exploration behaviour (black lines, mediolateral), and rolling behaviour (black lines, dorsoventral; modes ordered by ascending wavelength). B) Eigenvalues of the difference operators in our effective theory predict increasing frequency and damping, and decreasing variance, at shorter wavelengths. C) A large proportion of the variance in our datasets is explained by the longest wavelength modes (forward crawling, orange; backward crawling, green; unbiased exploration, black; rolling, purple). Equipartition of energy closely predicts the proportion of variance during crawling and unbiased exploration (blue). D and E) The frequency relationships between modes during crawling (D; top=forward, bottom=backward) and unbiased exploration (D, black) are well explained by our model (blue; 1-parameter nonlinear least squares fit to each behaviour; grey lines show experimental limits). F) Our effective theory closely predicts the average frequency relationship between the first axial and mediolateral modes during forward crawling and unbiased behaviour, respectively (blue lines), and correctly predicts that the average angular velocity of rolling should exactly match the average frequency of the first mediolateral bending mode during unbiased behaviour (p = 0.80, Mann-Whitney U-test, N = 36 bouts of rolling, N = 467 trials of unbiased exploration) but incorrectly predicts that frequency of forward and backward crawling should match (p = 0.009, Mann-Whitney U-test, N = 7 bouts of forward crawling, N = 5 bouts of backward crawling). We provide an example of how this forward/backward asymmetry can arise via self-interaction of the axial field in Appendix C.
There is an exact equivalence between modal analysis and PCA for mechanical systems in thermal equilibrium. Such systems obey the equipartition theorem which states that each mode has equal average energy 1 2β , where β = 1 k B T is the inverse thermodynamic temperature. This allows us to write the variance of the associated modal coordinate σ 2 i = Q 2 i in terms of the inverse temperature and the characteristic modal frequency.
Since the Hamiltonian contains no interactions between the modal coordinates, they should be uncorrelated when in thermal equilibrium (since for the modes to become correlated, they would have to interact in some way). This suggests that the covariance matrix for the mechanical system in its original non-modal coordinates can be written as where V is a matrix containing the mode shapes and Λ is a diagonal matrix containing the corresponding eigenvalues of the difference operators. Thus, for a mechanical system in thermal equilibrium the principal component vectors are exactly the mode shapes, while the variance of each principal component is inversely proportional to the corresponding eigenvalue. This tells us that, in thermal equilibrium, the proportion of variance explained by each principal component is entirely determined by equipartition, and depends only upon the eigenvalues of the difference operators, since we have In Appendices E, F, and G we exploit the structure of the difference operators to exactly analytically determine the mode shapes and their eigenvalues (which in turn determine the characteristic frequency scales, damping factors, and variance ratios according to Equations 8, 9, and 12 above) for the axial and torsional fields, while we approximate these quantities for the transverse fields using a matrix asymptotics approach.
First we consider the axial field. The circular symmetry introduced by imposing periodic boundary conditions on the axial field causes its modes to come in pairs with identical eigenvalues. The eigenvalues of the i'th pair of modes can be written as The value of the associated mode shapes (eigenvectors) at the k'th point in our spatial discretisation can be given as The zero eigenvalue must only be counted once in these expressions, since the corresponding "pair" of mode shapes are actually identical, and correspond to overall translation in the direction of the body axis. For the torsional field the modes are distinct, and their eigenvalues can be written as where there is again a single zero eigenvalue, corresponding to overall rotation around the body axis. The mode shapes (eigenvectors) are given as Each transverse field also has a set of distinct modes, although the mediolateral and dorsoventral fields have identical modes owing to their assumed circular symmetry. These fields have two zero eigenvalues, corresponding to overall translation and rotation in the coronal or sagittal planes, respectively. Our approximation for the i'th non-zero eigenvalue can be given as which becomes exact in the continuum limit N → ∞. The corresponding mode shapes are given by We compare the analytical modal analysis of our effective theory to PCA of real larval behaviour in Figure 2. We find that the proportion of variance explained by each mode is surprisingly well explained by equipartition of energy, with substantial variance explained by the first few long-wavelength modes. The first four axial modes explain 88.25% axial variance during forward crawling and 86.39% during backward crawling, compared to a theoretically predicted proportion of variance of 87.18% based on equipartition of energy. The first four transverse modes explain 97.97% of mediolateral variance during unbiased exploration and 99.84% during rolling, compared to a theoretical prediction of 96.83%. This suggests that shorter wavelength modes suffer from Boltzmann suppression -they have the same average energy as long-wavelength modes, but the same energy can produce much larger excitations in the long-wavelength modes than in their short-wavelength counterparts. The large amount of variance explained by the long-wavelength modes is also in agreement with the considerable spatial correlations we observed during larval behaviour.
At long wavelengths we obtain a strikingly good fit to the principal components of axial compression/expansion extracted from both forward and backward peristalsis, as well as to the mediolateral transverse principal components extracted from unbiased behaviour and the mediolateral-dorsoventral transverse principal components extracted from rolling behaviour. Since transverse bending was measured from a top-down view of the larva and rolling involves rotation of the larva relative to its substrate, the principal components extracted from rolling correspond to both mediolateral and dorsoventral bending (thus in Figure 2 we plot the principle components extracted from rolling behaviour alongside the "dorsoventral" modes and principal components extracted from unbiased behaviour alongside the "mediolateral" modes). Across all behaviours the fit is particularly close for the first two long-wavelength modes, which explain > 60% of axial variance during crawling, and > 94% mediolateral variance during unbiased behaviour. The first mode alone explains 96.78% mediolateral variance during rolling behaviour. One of the striking results of our analysis is that the eigenmaggot shapes are largely conserved between unbiased behaviour and rolling behaviour -our theory gives a strong explanation of why this must be the case (the shapes are encoded in the effective physics of the midline), but it would otherwise be highly non-intuitive. For instance, it may be surprising that a large amount of mediolateral variance (86.82%) during unbiased behaviour is explained by the characteristic C-bending mode that is most visually obvious during rolling, however our effective theory tells us that this must be true simply because equipartition of energy amongst the modes will always lead to a large amount of variance being explained by C-bending.
Moving on to consider the dynamics of larval behaviour, the relative frequencies determined by the eigenvalues of the difference operators provide a good fit to the average relative frequencies of all six axial characteristic frequencies observed during forward and backward crawling as well as the average relative frequencies of the first four mediolateral bending modes observed during unbiased behaviour. However, the absolute frequencies of forward and backward crawling in the real larva differ [5,37], requiring us to introduce different overall timescales for the two behaviours. We provide an example of how this forward/backward asymmetry may arise via an axial self-interaction in Appendix C.
We determine the average frequency ratio between the first mediolateral bending mode during unbiased behaviour and the first axial mode during forward crawling to be ≈ 0.52. In Appendix D we predict a ratio of 0.5 based on the inclusion of a prototypical nonlinear axial-transverse interaction term into our effective theory (based on the nonlinear curvature measure used in [30]), that leads to the transfer of energy between axial and transverse modes via parametric resonance.

Explanatory models of the larval behavioural repertoire
Motivated by the Boltzmann suppression of high-wavelength modes in our effective theory, we next aim to build low-dimensional explanatory theories of several larval behaviours by focusing only on the relevant long-wavelength excitations. In each case, we build a behaviour-specific theory with a minimal number of degrees of freedom corresponding to the first one or two relevant long-wavelength modes.
While our aim is not necessarily to explain the fine details of these behaviours, rather to provide a theoretical tool to aid in their conceptualisation and deeper understanding, we shall see that several of our "qualitative" theories actually provide an extremely good quantitative fit to the behaviour of the real animal.

Self-righting behaviour
We first turn our attention to the self-righting (SR) behaviour of the larva. This behaviour can be observed by experimentally preparing larvae in an "upside-down" configuration with their dorsal surface in contact with the substrate. Starting from this configuration, the larva appears to twist its body so as to attach its mouth hooks to the substrate, before propagating a wave of twisting along its body that brings the rest of the body fully into alignment with the substrate.
Our first observation is that the passive mechanical equilibrium of our effective theory, i.e. the minimum energy configuration, has all segments aligned with one another, since in this configuration there is no energy stored in the relative twisting or bending of body segments relative to one another. Thus, if one segment is able to attach itself to the substrate (for instance, via the mouth hooks at the head), and the larva relaxes all applied forces, all other segments will then passively align themselves and the larva will have achieved SR.
Since the symmetry proporties of our effective theory predict that torsional mechanics should be governed by the (damped) classical wave equation, the passive approach to the aligned equilibrium configuration should indeed occur via wave-like motion.
The primary problem facing the larva is thus how to align the head with the substrate, in order to achieve mouth hook attachment. Starting from the initial un-righted configuration, the larva will be able to align its head with the substrate by producing a torsional deformation along its length -this will rotate the mouth hooks around the midline, towards the substrate. A low-dimensional model of this deformation and the A low-dimensional model of self-righting behaviour. A, i) The larva begins in an "upside-down" equilibrium configuration. A, ii) muscular torques produce a distributed twist along the length of the body. A, iii) This allows the mouth hooks to attach to the substrate, producing a new stable equilibrium configuration. A, iv) Muscular torque is released, and the body passively un-twists towards its new equilibrium. A, v) The entire body is aligned with the substrate. B, top) Generalised torque applied to the first torsional mode. B, middle) First twisting mode -twist develops until mouth hook attachment, then relaxes. B, bottom) Average rotation of body segments relative to the substrate -no overall rotation occurs until after mouth hook attachment. C, top) Rotation of each segment relative to the substrate. All segments start "upside-down", then head and tail twist in opposite directions until mouth hook attachment, then all segments relax towards the aligned equilibrium configuration. C, bottom) Time taken for segments to align with the substrate, within some tolerance δ, reveals a travelling wave of alignment that propagates from head to tail. Blue verticle line indicates mouth hook attachment. Shown behaviour is valid for critically damped or overdamped torsion, ζ t,1 ≥ 1.
subsequent alignment can be constructed by retaining only two degrees of freedom from our effective theory -the overall rotation of the larva relative to the substrate, and the first torsional modal coordinate. This twisting mode characterises the amplitude of a monotonic, long-wavelength twist that is distributed along the entire length of the larva, and involves the head and tail twisting in opposite directions relative to one another. We can write the dynamics of the angle of the head relative to the substrate as a damped oscillator where ω t,1 , ζ t,1 are the characteristic frequency and damping factor of the first torsional mode, θ 0 is the initial non-righted orientation of the head (which we set to θ 0 = π to represent complete inversion of the body), and tau is a muscular twisting torque applied by the musculature. Provided the larva can sense its overall rotation relative to the substrate θ 0 , it can achieve head alignment by providing a constant twsting torque τ = −ω 2 t,1 θ 0 . Indeed, substituting this torque into Equation 19 gives so that the equilibrium configuration, for whichθ h =θ h = 0, clearly corresponds to θ h = 0, i.e. as the mechanics equilibrates to the applied torque, the head will align with the substrate. Following mouth hook attachment, the subsequent alignment of the body is governed by passive relaxation of the first torsional mode towards equilibrium, which is given by exactly the same dynamics as Equation 20 but with the modal coordinate Θ 1 substituted in place of θ h (see Appendix H).
Thus both the initial alignment of the head with the substrate and the subsequent alignment of the rest of the body with the head are governed by identical damped oscillator dynamics. The exact trajectory taken during SR in our model is thus governed by the characteristic frequency and damping factor of the first torsional mode (ω t,1 and ζ t.1 ). The role of these parameters in determining the behaviour of a damped oscillator is very well understood [55], and can be characterised as follows • Overdamped response (ζ t,1 > 1) -the system approaches equilibrium without any overshoot, with a double exponential timescale given by ζ t,1 ω t,1 ± Z 2 t,1 − 1ω t,1 . Increasing ζ t,1 or decreasing ω t,1 leads to a slower approach to equilibrium (i.e. slower SR), while decreasing ζ t,1 or increasing ω t,1 leads to a faster approach (faster SR).
• Underdamped response (ζ t,1 < 1) -the system approaches equilibrium and overshoots, followed by decaying oscillations ("ringing") around the equilibrium. The ringing frequency is 1 − ζ 2 t,1 ω t,1 , while the decay to equilibrium occurs at rate ζ t,1 ω t,1 . Increasing either ζ t,1 or ω t,1 leads to a faster approach to equilibrium (faster SR) while decreasing either parameter leads to a slower approach to equilibrium (slower SR).
• Critically damped response (ζ t,1 < 1) -the system approaches equilibrium in the least possible time, with no overshoot and a single timescale ω t,1 .
To summarise, our model of SR has two stages. First, a generalised force of constant magnitude is applied along the first torsional mode. This effectively shifts the equilibrium configuration and gives rise to a distributed twist along the body axis, causing the mouth hooks to align with the substrate. Second, the mouth hooks attach to the substrate and the larva releases the torsional force, so that the equilibrium configuration is shifted "back" again and consists of all segments coming into alignment. The release of the torsional force thus causes the larva to fully align with the substrate. The timescale for both stages of SR is dominated by the characteristic frequency and damping factor of the first torsional mode. Increasing the characteristic frequency (for instance, by increasing the shear modulus G of the larval cuticle or through neuromuscular effects, see Appendix B) should always lead to faster self-righting. Meanwhile, for a given characteristic frequency, self-righting will proceed fastest when the torsional mode is critically damped, for which the effective damping factor is equal to unity.

Rearing and hunching behaviour
Our analysis of self-righting showed that a simple constant input to the first torsional mode was sufficient to explain the low-dimensional dynamics of the behaviour. We note that similar constant inputs to the first dorsoventral bending mode or first axial mode produce deformations of the body that are qualitatively similar to rearing (provided the tail/posterior abdomen is attached to the substrate) or hunching, respectively ( Figure 4). In both these cases, and in our model of self-righting, the function of the constant input is to shift the equilibrium of our theory so that it corresponds to a new, deformed configuration. While we have not measured the kinematics of these behaviours we believe there is a strong likelihood that they will be well explained by low-dimensional dynamics in these modes; this provides a basic prediction for future experimental studies. Low-dimensional models of rearing and hunching behaviours. Starting from an initially undeformed configuration in which the tail/posterior abdomen is attached to the substrate (left) a constant-magnitude generalised force along the first dorsoventral bending mode will shift the passive equilibrium of the model to produce a "rearing" motion (top) while a constant-magnitude generalised force along the first axial compression mode will shift the equilibrium to produce a "hunching" motion (bottom).
Note that the roles of the effective mechanical parameters in our model of hunching and rearing are identical to their roles during self-righting, since all of the modes are governed by equivalent damped oscillator dynamics.

Peristalsis and rolling behaviour
Next, we turn our attention to peristalsis and rolling behaviours. In contrast to self-righting, which has a clearly defined "end-point" in which the body is in equilibrium, these behaviours involve a continuous, ongoing, seemingly out-of-equilibrium, motion of the body.
For instance, in peristalsis, a wave of segmental compression and expansion passes along the body. In forward crawling, this wave travels from the tail of the animal, along its anteroposterior axis, until it reaches the head. At this point the head, tail, and internal viscera move in synchrony [5], and a new wave is initiated at the tail. In backward crawling the behaviour is similar, except the wave travels from head to tail and propagates at around half the speed of its forward-propagating counterparts ( Figure 2). In both cases, the wave has a single "peak" of maximal compression, but is spread across several segments, so that it appears largely as a broad sinusoid of compression and expansion (recall that the first two axial modes explain around 60% of variance during these behaviours). In rolling, the larva enters a C-bend configuration (which looks a lot like the first bending mode of our model -recall that this mode explains around 95% of variance during the behaviour). The C-bend stays in the plane of the substrate while the body rotates, so that the C-bend propagates around the body axis in a clockwise or counterclockwise direction. For clockwise rolling, the C-bend might first be directed to the right of the body, then ventrally, then to the left, then dorsally, then back to the right; for counterclockwise rolling the sequence is reversed.
To the eye, rolling and peristalsis behaviours could scarcely look more different. However, we note that we have already seen some similarities between the two behaviours. First, both behaviours involve longwavelength modal excitations, with the first two modes explaining substantial variance. Crucially, in both behaviours, these long-wavelength excitations have degenerate (identical) frequencies, owing to the underlying circular symmetries of the axial and transverse fields that was introduced during the construction of our effective theory -this degeneracy is clear in Figure 2, where there is a "pairing up" of frequenciesthe first two axial modes are identical, as are the frequencies of the first mediolateral and first dorsoventral bending mode, and the pairing continues for the higher modes. This degeneracy is the key to our analysis of rolling and peristalsis, as it enables us to study both behaviours using an identical Hamiltonian formalism.
In particular, we take our Hamiltonian expressed in modal coordinates (Equation 7), keep only the terms for a pair of modes with identical frequency, and transform to polar coordinates (write Q 1 = A cos φ for the first mode and Q 2 = A sin φ for the second, and then transform the momenta; see Appendix I for more details).
We are left with where ω is the frequency of the two modes, A is the "amplitude" of the modes, given by A = Q 2 1 + Q 2 2 , p A is the canonical momentum associated with this amplitude, and M is the canonical momentum associated with the "phase" angle, given by φ = tan −1 Q 2 /Q 1 . Note that the phase φ does not turn up in the Hamiltonian at all. This represents the underlying circular symmetries that were introduced during the construction of our theory -the Hamiltonian is symmetric under rotations of the configuration space. To see this, we note that such a rotation amounts to adding a constant to the angle φ, but since φ is absent, the Hamiltonian does not change. A beautiful result from classical mechanics known as Noether's theorem tells us that this symmetry must be tied to a conserved quantity. Indeed, Hamilton's equations (Equation 6) tell us that (in the absence of dissipation)Ṁ = ∂H ∂φ = 0 (22) so that the momentum M is conserved -it stays constant over time.
Let us pause to consider the biological meaning of these observations. Imagine a situation in which the larva is in a C-bend configuration as in Figure 5. In our effective theory, the energy associated with the bend does (right) our statistical mechanics model closely predicts the distributions of larval kinematic variables for the first two axial modes observed during forward crawling (blue). The single free parameter for these distributions, inverse thermodynamic temperature β, was estimated by fitting the amplitude distribution using nonlinear least squares.
not depend on how it is oriented relative to the body in a dorsoventral-mediolateral coordinate system -a dorsal bend has the same energy as a lateral bend of the same amplitude, and bends at intermediate angles also have the same energy. These energies are equivalent precisely because the characteristic frequencies of dorsoventral and mediolateral C-bending are identical. We call this a rotational symmetry because the numerical value of the energy does not change as we rotate the bend around the body. Noether's theorem then tells us that there is a conserved quantity associated with this symmetry -in our case it is the angular momentum associated with propagating the bend around the body. In the case of the axial modes, the "rotation" corresponds to moving a pattern of compression/expansion along the body and "through" the visceral piston that mechanically links the head and tail -this is a symmetry of our effective theory because the energy of the compression/expansion will be the same regardless of how it is shifted through the body. Noether's theorem then tells us that a generalised form of angular momentum associated with "rotating" this compression through the body is conserved. We illustrate both of these examples in Figure 5.
Since the momentum M does not change, we can treat it as a fixed parameter in the Hamiltonian (Equation 21), so that we are left with the motion of the amplitude coordinate A in an effective potential energy U ef f that depends upon the particular, constant, value of M (this reduction in dimensionality by introducing an effective potential is a common trick when solving problems with a rotational symmetry [55]). When M is non-zero, the effective potential has a minimum at A = M/ω, signalling that the new minimum-energy Metropolis samples=dark blue; analytical distribution=light blue, real distribution=black). The single free parameter for these distributions, inverse thermodynamic temperature β, was estimated by fitting the amplitude distribution using nonlinear least squares. equilibrium solution, taking into account conservation of the momentum M , has a non-zero expectation value -the equilibrium solution now corresponds to continuous propagation of sinusoidal compression waves through the body or C-bends around the body axis. Indeed, the equilibrium solution has A = M/ω so that the Hamilton's equations tell usφ This means that the phase angle evolves as φ = ωt + C, and the minimum-energy solution has the modal coordinates evolving sinusoidally at the characteristic frequency, This provides the a novel experimental test of our theory of rolling. Interpreting Q 1 , Q 2 as measuring the extent of mediolateral and dorsoventral C-bending, respectively, this equation tells us that the C-bend must complete one rotation around the body axis in a time period 2π/ω b,1 where ω b,1 is the characteristic frequency of C-bending. One rotation of the C-bend about the body axis corresponds to one complete rotation of the larva about its body axis, since the C-bend stays fixed in the plane of the substrate. Therefore, our effective theory makes the surprising prediction that the mean angular velocity of rolling should be equal to the mean timescale of the first mediolateral bending mode during unbiased behaviour. Indeed, we have measured the mean angular velocity during n = 36 trials of rolling behaviour and measured the mean timescale of the first mediolateral bending mode during n = 467 trials of unbiased behaviour, and find no statistically significant difference in the measurements (p = 0.80, Mann-Whitney U-test).
Armed with our low-dimensional theory of peristalsis and rolling, and motivated by the fact that an appeal to thermal equilibrium was able to correctly predict the principal components, and their proportions of variance, during real larval behaviours, we next ask whether equilibrium statistical mechanics is able to predict the distributions of larval kinematic variables during rolling and peristalsis. For calculational simplicity, we use the canonical ensemble of equilibrium statistical mechanics, in which we constrain the average energy to be constant and treat the momentum M as an additional constraint. In this case the probability of a given mechanical state is given by a modified Boltzmann distribution (see Appendix J), which can be written as follows where H(state) is the energy of the mechanical state (given by the Hamiltonian, Eqn 21), M (state) is the state's degenerate angular momentum, which is constrained to take a fixed value M 0 due to the Dirac delta function δ, β is the inverse thermodynamic temperature, and Z is the partition function which can be determined by normalising the probability distribution. We have calculated the partition function in Appendix J, and can write it as follows where is the reduced Planck's constant. Note that the appearance of the absolute value of the momentum |M | in the partition function signals non-analyticity of the free energy; we explore the consequences of this non-analyticity in the following section. For our system the Hamiltonian can be written as a sum of a kinetic and (effective) potential term, so the Boltzmann distribution factors, allowing us to easily compute the marginal distribution over configuration space in terms of the effective potential. The configuration space distribution can be written as where Z c is the configurational partition function, which can again be determined by normalisation. The predicted distributions of larval kinematic measurements can then be predicted directly from this distribution. Note that the minimum-energy solution with A = M/ω now has the highest probability, since the Boltzmann distribution weights low-energy states more highly, so that the minimum-energy equilibrium solution corresponds to the expected or "average" behaviour.
In Figures 6A and 7A we plot the low-dimensional minimum-energy behavioural trajectories predicted by our model (orange) compared to real measurements (black) for forward crawling and rolling behaviours, respectively. Figure 6B and 7B show the trajectories in the modal configuration space -here the minimumenergy trajectory (orange) forms a perfect circle of radius M/ω, while the real data (black) and the equilibrium statistical mechanics model (blue) show a spread around this trajectory. The predicted and measured marginal distributions over the individual modal coordinates, the phase angle, and the amplitude are shown in Figures 6C and 7C. Note that for rolling behaviour we have only been able to experimentally measure the amplitude, and not the other kinematic variables -our theory thus produces novel predictions for future experimental investigations of rolling.

Unbiased substrate exploration behaviour
We next return to our expression for the partition function of a pair of modes with identical frequency, in thermal equilibrium at an inverse thermodynamic temperature β and with conserved momentum M (Equation 26). As we noted earlier, the partition function depends upon the absolute value of M and so is non-analytic, suggesting that the Helmholtz free energy should also be non-analytic; such non-analyticities very generally represent the presence of phase transitions [49]. The Helmholtz free energy is given by which is, as expected, a non-analytic function of M . Indeed, the first derivative of F with respect to M , which gives the mean frequency of rolling or crawling, is discontinuous at M = 0, since This indicates the presence of a first-order, abrupt phase transition as M goes through zero [49]. For M < 0 we have counter-clockwise rolling or backward crawling and for M > 0 we have clockwise rolling or forward crawling. The discontinuity of the free energy indicates that the forward/clockwise and backward/counterclockwise behaviours cannot co-exist within our effective theory -they are separated by a sharp phase transition.
Let us consider what happens exactly at the transition (M = 0), focusing on the transverse modes. When the degenerate angular momentum vanishes, the transverse bending must be confined to a plane, as the bends cannot rotationally propagate around the midline (if they did, we would have M = 0). This is essentially what happens during unbiased substrate exploration behaviour, during which the larva primarily moves parallel to the plane of the substrate [30,63,64]. We therefore asked whether our Boltzmann distribution could predict the statistics of unbiased exploration in the case M = 0. Assuming that all bending occurs along the mediolateral direction, we can write the Boltzmann distribution for the i'th mediolateral modal coordinate Y i as with partition function Thus the predicted modal distributions are simply zero-centred Gaussians with variances 1/βω 2 b λ i , in line with our earlier results using the equipartition theorem. Since the eigenvalues λ i are fixed properties of the difference operators in our discretised effective theory, and we have already fit the parameter ω b to match the characteristic frequencies observed during larval behaviour, we are left with a single free parameter β with which to fit all nine observed modal distributions. Figure 8 shows the results of this fitting procedure. There is an excellent agreement between theory and experiment, with our model closely predicting the variances of all 9 distributions ( Figure 8C), and closely matching the actual form of the distributions ( Figure 8B), although the observed distributions have slightly heavier tails than our predicted Gaussian distributions.
Since our mediolateral modal coordinates form a complete basis with which to describe the planar postures of the larva, we are also able to use our distributions (Eqn 30) to predict any other measurement depending only on the posture of the animal. To illustrate this, we compute the distribution of 3-point body bends for the larva, a common experimental measure in the larval literature (the 3-point body bend is the angle formed between a line joining the tail and midpoint and a line joining the midpoint and head). There is again a striking agreement between the prediction of our effective theory and the measured distribution, up to very large bending angles ( Figure 8D).

Discussion
In this paper we have presented an idealised model of larval midline mechanics. It is the first model to describe the fully 3-dimensional motions of the larval midline, including stretching along, bending of, and twisting about the body axis. Owing to the simplicity of the ideal model, we are able to derive many results analytically without having to resort to simulation. Indeed, the model makes strong quantitative predictions for the dynamics of peristalsis, rolling, and self-righting behaviours, and also predicts statistical properties of these behaviours and of unbiased exploration behaviour.
We constructed our long-wavelength effective theory from first principles, using a Lagrangian expansion approach borrowed from modern physics. We believe that this approach is very powerful -it demonstrates the close ties between symmetry properties of the larval body and the dynamics and statistics of the animal's behaviour. To our surprise, symmetry considerations combined with basic tools from classical mechanics and equilibrium statistical mechanics are sufficient to explain a wide range of the animal's behaviour. This is because our imposed symmetries fix the difference operators in our theory, and these in turn determine the dominant patterns (called modes, principal components, or eigenmaggots) in larval behaviour, along with the relationships between their characteristic timescales. Furthermore, the relative dominance of each pattern (the proportion of variance associated with each eigenmaggot) appears to be predicted by simple equipartition of energy. This means that all eigenmaggots have the same energy, on average, leading the longwavelength eigenmaggots to dominate larval behaviour via Boltzmann suppression; the same average energy can produce a larger average excitation of the long-wavelength eigenmaggots than the short-wavelength eigenmaggots.
The statistical suppression of short-wavelength modes allows us to build low-dimensional descriptions of larval behaviours, by focusing only on the long-wavelength modes (the first few eigenmaggots). Since each mode has deterministic dynamics governed by the extremely well understood damped harmonic oscillator equation, we are able to build remarkably simple models of self-righting, rearing, and hunching, in which constant modal force inputs shift the midline's passive equilibria. Our low-dimensional analysis of rolling and peristalsis is based on the O(2) circular rotational symmetries of the axial and transverse modes in our effective theory; this shared symmetry property makes the physics of the two behaviours surprisingly similar.
To analyse the behaviours we use equilibrium statistical mechanics -this allows us to easily study the effects of a give-and-take of energy into the body via the musculature without needing to specify the exact form of the neuromuscular dynamics, and is motivated by our earlier success in applying equipartition to predict the proportion of variance in each mode. The Boltzmann distribution then provides an excellent fit to real larval statistics during rolling and peristalsis, provided we take account of the additional conserved momentum signalled by the O(2) rotational symmetry. From here, the averaged peristalsis and rolling behaviour of the larva is shown to correspond to a classical, sinusoidal, conservative trajectory, at the minimum of an effective potential energy.This deterministic trajectory provides an excellent fit to the real low-dimensional trajectories during rolling and peristalsis.
We are finally able to show that the statistics of unbiased exploration behaviour appear at an abrupt phase transition between clockwise and counterclockwise rolling. During rolling, bends propagate around the axis of the body, associated with a certain momentum M . The phase transition occurs for M = 0, when the bends no longer propagate (as in rolling) and are instead confined to a plane (as in unbiased behaviour). At the transition, the predicted behavioural statistics go through a qualitative change -the predicted distributions over the bending modes switch from being non-Gaussian, non-zero-centred, with a positive skew, to being zero-centred Gaussians. Amazingly, the new statistics that emerge at the transition between clockwise and counterclockwise rolling very closely match the observed statistics of unbiased behaviour, with a very close fit to all nine observed bending mode distributions. We are able to use this statistical description to predict the distribution of 3-point body bends during unbiased behaviour (a common measure used in the larval literature [65]), up to very large body bends.
It may seem odd that our theory is able to account for the effects of the neuromuscular system without describing any details about how the nervous sytem is actually instantiated in the larva, and by only describing the physics of the midline. We can build some intuition for why our theory is able to accomodate neuromuscular effects by noting that the behaviours of the animal consist of the motion of a physical object, the animal's body, and that the role of the neuromuscular system in producing behaviour is to exert forces on this object. Crucially, the effects of muscular forces do not depend upon how they were produced -they do not depend on the neural mechanisms used to coordinate force production -they only depend upon the relationship between the force and the mechanical state of the object upon which they act. In particular, depending upon the relationship between the force and the mechanical state, muscles may function in a variety of different effective mechanical roles. [1] summarises this by saying that muscles may act as motors, brakes, springs, or struts. The distinction between these roles is essentially whether the musculature provides power flow into the body (muscles act as motors), power flow out of the body (muscles act as brakes), or provides power flow into and out of the body with no average change in the mechanical energy (muscles act as springs or struts). Thus it should not be too surprising that neuromuscular effects can be described within a purely physical theory, either as modifications to the effective Rayleigh dissipation function (power flow in or out of the body, i.e. the effects of muscles as motors and brakes) or the effective Lagrangian (no average power flow, i.e. the effects of muscles as springs or struts). In Appendix B we provide a more detailed account of the incorporation of neuromuscular forces, showing that in the two particularly important cases of forces that explicitly depend on the mechanical state (e.g. reflexes) or periodic forces that do not explicitly depend on the mechanical state (e.g. feedforward control via CPGs without sensory feedback), the overall effect at long length and time scales can be almost entirely accomodated by a simple rescaling of the phenomenological parameters.
The fact that equilibrium statistical mechanics can explain observed probability distributions for larval kinematics measurements is genuinely surprising. We note that the application of equilibrium statistical mechanics to our effective theory yields a model that is very similar to a discretised form of the Worm-Like Chain model that is commonly used in polymer physics [53,54]. Thus, our results suggest that the statistics of macroscopic bending of the larva during unbiased exploration are similar to those of the microscopic thermal "jiggling" of a polymer chain.
Crucially, the application of statistical mechanics to our mechanical model necessarily presupposes some mechanism that is not contained directly in our effective theory. This should be clear since the effective theory gives rise to linear deterministic dynamics, and so contains no mechanism by which trajectories may "wander" or "mix" in the phase space -in other words, the deterministic dynamics of the theory cannot explain the statistics of the theory without modification. This is actually a fairly fundamental issue in statistical mechanics, not specific to our theory. It is usually overcome by appealing to one of several (related) mechanisms [49,66]: • microscopic chaos, in which the microscopic laws of a theory are simply presupposed to give rise to sufficiently ergodic chaotic exploration of energy surfaces, permitting us to treat all points on an energy surface as having a uniform probability; surprisingly this often gives correct results even when the underlying physics is provably non-chaotic.
• weak coupling in which the separate subsystems in a theory (corresponding to modes/eigenmaggots, in our case) are able to transfer small amounts of energy between themselves.
• introducing a heat bath which allows transfer of energy between the system and its "environment" (which might include some weakly coupled subsystems, and might also capture the effects of the neuromuscular system, as described below).
It is quite possible that small modifications to our effective theory may lead to chaotic dynamics or energetic coupling. For instance, our previous work on the planar mechanics of unbiased exploration [30] suggests that beyond the low energy limit we have considered here, it is very likely that the mechanics of bending and stretching of the midline should become energetically coupled, leading to chaotic motion. However, it is currently unclear whether this model is sufficiently ergodic to explain the appearance of equilibrium statistics, and in particular whether the model reproduces equipartition of energy. Recent work in the nematode C elegans has experimentally demonstrated the presence of a symmetric Lyapunov spectrum for the animal's dynamics, strongly suggestive of deterministic chaos arising due to (damped, driven) Hamiltonian mechanics; if such a spectrum is present in larval behaviour, it may go some way to explaining the animals' behavioural statistics.
The success of equilibrium statistical mechanics and Lagrangian field theory in predicting the statistics of larval behaviours suggests that the action of the neuromuscular system can be approximated by three effects: (1) modifications to the dissipative Lagrangian effective field theory, for instance due to reflexes or CPG entrainment introducing new forces that effectively depend upon the fields (modelled by effective energetic or dissipative terms), (2) maintenance of conserved quantities of the effective theory (for instance maintaining the rolling angular momentum in our theory), and (3) introduction of a "heat bath" to add non-specific energy to the system. In other words, the effect of the neuromuscular system on the scale of aggregate behaviour appears to be adequately explained by a statistical mechanics model that makes absolutely no reference to the particular details of the nervous system. However, this clearly does not mean that our effective theory has nothing to say about the neuroscience of larval behaviour. Indeed, we believe that our theory may be able to provide mechanistic hints for experimental biologists. As an example, consider a high-throughput experiment in which a certain effect on the aggregate statistics of mutant larvae can be explained as a modification of the effective stiffness of the animal, it suggests that the mutation may effect either elastic properties of the cuticle or the entrained neuromechanical phase relationships (since forcing inphase with position alters the effective stiffness, as muscles begin to act as springs [1]). Thus, by abstracting away from any neuroscientific details, our theory is able to provide an excellent explanation of the overall effects of the nervous system on behaviour, and this explanation may in turn have mechanistic relevance.
Given the simplicity of our starting assumptions we believe our effective theory may be readily applicable to studying the behaviour of other animals with similar body plans to the larva. This includes several other important "model species" such the nematode worm C. elegans, the only animal for which a complete neural connectome is currently available, and the zebrafish, thus forging a stronger link between invertebrate and vertebrate systems.

Author roles and acknowledgements
Jane Loveless developed the theoretical aspects of this work, conducted all data analysis, and wrote the bulk of the manuscript. Alastair Garner conducted experiments on rolling and unbiased behaviour. Raouf Issa conducted experiments on self-righting. Barbara Webb provided guidance on the development of the theory and aided with editing the manuscript. Tomoko Ohyama supervised collection of experimental data on rolling and unbiased behaviour, conducted experiments on rolling and unbiased behaviour, and aided with editing the manuscript. Claudio Alonso supervised collection of experimental data on self-righting behaviour, provided guidance on writing the manuscript, and aided with editing the manuscript. Tomoko Ohyama and Claudio Alonso together provided monetary support for Jane Loveless during the research underlying this paper.
Data for forward/backward peristalsis was provided by Stefan Pulver. This data was previously publised in [37].
Illustration of self-righting behaviour in Figure 3 was produced by Eleanor Savage.
Jane Loveless would like to thank Kostantinos Lagogiannis, Greg Stephens, Tosif Ahamed, and Antonio Carlos Costa for fruitful discussions during the development of this paper.

A Detailed construction of the effective theory
To construct our effective theory of the larval AP-axis (midline), we will essentially perform a search in "theory space", by first writing the most general possible Lagrangian and Rayleigh densities and then eliminating any parts of them that do not meet certain physical requirements. This is similar to a variety of approaches taken in many branches of modern physics [3,49,50,58].
Our key requirements during the construction of our low-energy effective theory are analyticity, Galilean symmetry, and reflection symmetry of the Lagrangian. We shall also require that our theory avoids the Ostragradsky instability, which forbids the inclusion of time-derivatives of higher than first order in the Lagrangian, and that our Lagrangian density includes a "standard" kinetic energy term proportional to the sum of squares of the field velocities. This last condition rules out a single term, corresponding to the Rayleigh correction that would account for the rotational inertia associated with the larva's cross-section; this correction should predominantly effect the high-frequency behaviour of the larva so it is neglected. We will see that when we focus our attention to long length and time scales (equivalently low energy scales) we are left with only one possible effective Lagrangian field theory of the midline.
We begin by introducing a set of four scalar fields x(s, t), y(s, t), z(s, t), and θ(s, t), parametrised by time t and the arc length s along the undeformed axis of the body, and respectively measuring the axial, mediolateral, dorsoventral, and torsional displacements of the midline away from an undeformed configuration of the body (Figure 1).
We then write a completely general Lagrangian density for the midline in terms of the scalar fields and their derivatives with respect to time and space.  Table 1. All candidate terms in the fourth-order single-field Lagrangian density expansion. Each candidate term is given in terms of a general field φ, along with the order of the derivative (δ), amplitude scale ( ), and characteristic time and space scales (τ and σ, respectively). We then list whether each term satisfies the requirements of our effective theory (x indicates satisfaction). The axial and torsional Lagrangian densities must be Ostragradsky stable (OST), invariant under field translations φ → φ + ∆ (TRAN), and invariant under reflections φ → −φ (REF). We also require that the transverse Lagrangian densities have a standard kinetic energy term (∂ t φ) 2 (KIN), which rules out the Rayleigh correction term (∂ ts φ) 2 and the relativistic correction (∂ t φ) 4 , and they must be invariant under infinitesimal field rotations φ → φ + ∆s (ROT). Intuitively, the axial and torsional Lagrangian densities can only include terms which have x's in all three columns OST, TRAN, REF and the transverse Lagrangian densities can only include terms which have x's in all five columns OST, TRAN, REF, KIN, ROT. From the remaining terms for each field, we include only those of lowest order in the space and time scales (lowest τ + σ). For instance, for the axial field x we keep only the second-order terms (∂ t x) 2 and (∂ s x) 2 , since they will dominate the (otherwise allowed) fourth-order terms (∂ t x) 4 , (∂ s x) 4 , (∂ ts x) 2 , and (∂ ss x) 2 on long length and time scales.
For now we will focus on the separate densities for each field; we will return to consider the interaction Lagrangian density L I towards the end of this section. For concreteness, we will consider the axial Lagrangian density and we clearly have similar expressions for the other Lagrangian densities.
We then assume that each Lagrangian density can be written as an analytic function of its associated field and the derivatives of this field. Focusing specifically on the axial Lagrangian density, this means we can perform the expansion Where f (x) and the factors g(x), h(x), etc are arbitrary functions of the axial field. We characterise each term in the expansion according to its "order". To illustrate this procedure, we consider the term, where we have discarded the arbitary field-dependent factor (we will see shortly that these factors are unimportant). We denote the order of the derivative by δ, so for this term we have δ = 4 since it contains a fourth derivative. We introduce an amplitude scale by multiplying the field by a constant, φ → a · φ, and introduce time and space scales by similarly multiplying the time and arc-length by constants t → t · t 0 , s → s · s 0 . Substituting these scaled quantities into the term above, we arrive at We then denote the exponent of the amplitude scale a by , so for this term we have = 2 and we say that it is of second order in the amplitude. We similarly denote the exponent of t 0 by τ and the exponent of s 0 by σ. Thus for this term we have τ = 4 and σ = 6, and we say that the term is of fourth order in the time scale and sixth order in the spatial scale. On long length and time scales, and small amplitude scales, terms characterised by lower values of , τ , and σ must dominate the Lagrangian density. In Table 1 we provide an exhaustive list of all terms in the Lagrangian density expansion up to fourth order in δ, , τ , and σ.
Having explicitly written out our general Lagrangian density, we then begin to use our physical requirements to successively remove terms. Table 1 shows exactly which terms are ruled out by which requirements, up to fourth order.
First, we can remove all terms involving time derivatives of higher than first order from the kinetic energy densities, as the inclusion of such terms leads to the Ostragradsky instability [48].
We then consider the symmetry under overall translations and rotations, which represents our belief that the physics of the larva should not change depending upon its absolute position and orientation in space.

Translational symmetry implies that the Lagrangian density must be invariant under shifts in the fields
x → x + ∆, or y → y + ∆, or z → z + ∆. For this to hold generally, the arbitrary functions f , g, etc in the axial and transverse Lagrangian densities must be constants, independent of the fields. Since addition of a constant to the Lagrangian will not change the equations of motion [55], we can thus also safely remove the constant f and any term containing a power of the field alone, e.g. φ, φ 2 , φ 3 , φ 4 in Table 1.
We then look for terms obeying rotational symmetry. For the torsional field θ this simply requires invariance under field translations θ → θ + ∆.
To study rotations of the axial field x and transverse field y we start with the larva in its undeformed configuration, aligned with the X axis of an inertial Cartesian coordinate system, and with the origin of the X axis coinciding with x(0), the location of the tail. In this configuration we have x(s) = s and y(s) = y. A rotation through angle ∆ about the tail in the X, Y plane of the Cartesian system will transform the fields to new values x , y given by x (s) = s cos ∆, y (s) = y + s sin ∆ (38) For infinitesimal rotations ∆ → 0 this becomes x (s) = s = x(s), y (s) = y + ∆s = y(s) + ∆s (39) so that the axial field remains unchanged while the transverse field transforms as y → y + ∆s. We can clearly repeat our argument for rotations in the X, Z plane, so that the other transverse field must similarly transform as z → z + ∆s. Requiring invariance of the Lagrangian densities for the y and z fields under these transformations eliminates all terms involving spatial derivatives of less than second order (δ < 2), since their inclusion would lead to the appearance of terms proportional to ∆ in the Lagrangian densities under rotation (Table 1). Now we consider reflection symmetry. This reflects our belief that the physics of similar left-handed and right-handed bends should be equivalent, as should be the physics of similar dorsal or ventral bends, and similar clockwise or counterclockwise twists along the AP-axis. Reflection symmetry of the axial field simply represents that a motion in the anterior or posterior direction should be the same regardless of the overall direction the body axis is pointing in. It can be thought of as a special case of rotational symmetry combined with mediolateral or dorsoventral reflection -for instance we may perform an axial reflection by first rotating the body through 180 degrees in the mediolateral plane before performing a left/right reflection. Formally, reflection symmetry requires each Lagrangian density to be invariant under a discrete transformation of the fields as x → −x, y → −y, etc. This symmetry clearly rules out all terms raised to an odd power in the Lagrangian densities.
Our last requirement is that the Lagrangian densities should contain a "standard" kinetic energy term consisting of the square of the first time derivative of the field. This removes a term (∂ t φ) 4 , representing a relativistic correction to the kinetic energy that is unimportant at the velocities experienced by the larva (which are very small compared with the speed of light), and additionally eliminates a term (∂ ts φ) 2 , which corresponds to a correction to the kinetic energy of a beam first proposed by Rayleigh. Without the Rayleigh correction, the kinetic energy is composed only of the translational kinetic energy of the beam, with the mass considered as being concentrated exactly at the midline. The Rayleigh correction incorporates additional kinetic energy due to the distribution of mass across the cross-sections of the beam, and the rotation of these cross-sections as the beam bends. The primary effect of this correction is a lowering of the higher modal frequencies; since in our theory, and in the real larva, the majority of variance is explained by the lower frequency modes, we do not lose much by eliminating this term.
Finally, we restrict attention to long length and time scales, by retaining only those terms in the Lagrangian densities of lowest order in τ and σ. We are able to construct axial and torsional Lagrangian densities that are second order in both τ and σ, finding and The transverse Lagrangian densities cannot contain any terms involving the first space derivatives, due to the requirement of rotational invariance (see above and Table 1). The only term we can include is the square of the second space derivative of the field, which is fourth order in σ, so that the Lagrangian densities are given by and We now consider the interaction Lagrangian density L I , including terms up to second order in the time scale and up to fourth order in the length scale, since terms higher than this will be dominated by the terms we have already included in the Lagrangian densities for the individual fields. Clearly this interaction Lagrangian must simultaneously satisfy all of the requirements we placed on the separate Lagrangian densities for each scalar field. We note that the interaction Lagrangian cannot include any "self-interaction" terms involving only a single field and its derivatives, as such terms have already been accounted for during our previous discussion. Thus the interaction Lagrangian must consist of a sum over products of different fields and their derivatives only, where each term may be raised to an integer power.
The requirement of translation invariance immediately rules out the appearance of any term including the fields themselves, or powers of the fields, as a factor, such as xy, xz 2 θ, or θ ∂z ∂s , since such terms clearly are not invariant under general translations of the fields.
Reflection invariance further removes all terms which contain a factor which is not raised to an even power, since reflection of a single field in the interaction would cause these terms to change sign. Taken together, these requirements tell us that the lowest order interaction terms must consist of products of first derivatives raised to a power of 2. There are three such classes of interaction, the potential interaction where φ 1 = φ 2 denote the two interacting fields. Clearly the kinetic interaction is of fourth order in τ , so this interaction can be discarded from our Lagrangian density. Meanwhile, our requirement that the Lagrangian density have a "standard" kinetic energy term rules out any mixed interaction from appearing in the Lagrangian. Rotational invariance tells us that neither the mediolateral nor dorsoventral fields y, z may take part in a potential interaction since the Lagrangian cannot contain their first spatial derivatives. This means that our only viable interaction term is a potential interaction between the axial and torsional fields. However, we note that the axial-torsional interaction potential is of fourth order in σ while the separate axial and torsional potential energies are of second order in σ. Thus in the limit of low long wavelengths, we expect the axial-torsional dynamics to be dominated by their individual potentials, with a negligible contribution from their interaction potential. We therefore discard this interaction from the Lagrangian.
Since all other possible interaction terms are necessarily of higher order than the three we have just considered, and they are already of too high an order to include in the Lagrangian density, we conclude that in the limit of long length and time scales the four scalar fields must be effectively non-interacting.
We now derive the equations of motion for each field using the Euler-Lagrange equations, finding for the axial and torsional fields These are clearly recognisable as linear wave equations. For the transverse fields we have These are also clearly recognisable as Euler-Bernoulli beam equations. Together, our four equations of motion give the dynamics of linear elastic beam theory [57].
We next take finite difference approximations to the spatial derivatives in the equations of motion. To do so, we must first set boundary conditions. We set free-free boundary conditions on the transverse and torsional Lagrangian densities, and set periodic boundary conditions for the axial Lagrangian density, representing the coupling of axial head and tail motions via the visceral piston of the larva [18]. The discretised dynamics thus becomeẍ where the vectors x, y, z, θ contain the components of the scalar fields sampled at discrete points along the midline. The vectors y, z, θ each contain N components, while the presence of periodic boundary conditions on the field x reduces the number of components to N − 1, since head and tail displacements are equal. We have absorbed additional constants associated with the discretisation of space into the constants c x , c y , c z , c θ since these were essentially arbitrary. D 2,c is the (N − 1) × (N − 1) circulant second difference matrix (i.e. second difference matrix with periodic boundary conditions) D 2,f is the N × N second difference matrix with free-free boundary conditions and D 4 is the N × N fourth difference matrix with free-free boundary conditions These equations of motion can be derived from the corresponding Lagrangian as can be checked by the interested reader. We then construct the low-energy effective Hamiltonian via a Legendre transformation which is the equation that appears in the main text.
The low-energy effective Hamiltonian that we have derived can describe only the energy-conservative part of the midline mechanics. To account for the flow of energy through the body we now derive an effective Rayleigh dissipation function. Our Hamiltonian gives the total energy associated with a given state of the midline at an instant in time, while our Rayleigh dissipation function is related to the power flow associated with that state, at that instant. We will be able to derive the equations of motion for each discretised field from the dissipative Hamilton's equationsẋ where R is the Rayleigh dissipation function.
We proceed as before, by writing out the most general possible Rayleigh dissipation density for a continuous theory of midline mechanics, and then eliminating terms according to physical requirements. We will then limit ourselves to considering only those contributions that dominate at long wavelengths and long timescales, and will finally discretise in space. The key requirements will be Galilean symmetry (incorporating translational symmetry, rotational symmetry, and symmetry under boosts), reflection symmetry, and analyticity of the Rayleigh function. Our general Rayleigh dissipation density may be written where the first four terms are Rayleigh dissipation densities for the individual scalar fields while the final term allows for energy dissipation occuring due to the interaction of the fields. As before, we first focus only on the dissipation functions for the individual fields. Assuming analyticity, we have for the the axial dissipation function ∂ 2 x ∂t∂s , + . . . (63) As before, translation symmetry further tells us that the functions f , g, etc must in general be constants.
We can immediately rule out the occurrence of any terms depending only upon the fields themselves or their spatial derivatives, and rules out any terms involving higher time derivatives of the fields beyond first order. This is because the Rayleigh function only enters the dissipative Hamilton's equations via it's derivative with respect to the first time derivative of the field.
We can also rule out any term involving only the first time derivative of the field since they are not invariant under Galilean boosts of the form x → x + vt where v is the uniform velocity of the boost. This leaves us with only mixed space-time derivatives of first order in time. Our argument so far is equivalent to saying that the midline cannot dissipate energy simply due to its motion relative to our coordinate system, and that dissipation can occur only due to the relative motion of different parts of the midline.
We now apply reflection symmetry x → −x, y → −y, etc which again rules out all terms raised to an odd power. The Rayleigh function must therefore be of the form where the first term is second order in σ, the second and third terms are fourth order, and the final term is of eighth order. Restricting ourselves to long wavelengths, we keep only the second order term in the axial and torsional dissipation functions. For the transverse dissipation functions this term is not viable, as it is not invariant under (infinitesimal) continuous rotations of the form y → ∆sdt, z → ∆zdt; these rotational motions should not cause dissipation of energy in the absence of substrate interaction. The second term is also ruled out by this invariance requirement, so we keep only the third term, which is of fourth order in σ.
We note that, as for our interaction Lagrangian, the interaction terms in the Rayleigh dissipation function are once again necessarily of a higher order than the dissipation functions for the individual fields, so that we may rule out dissipative interactions between the fields at long length and time scales.
Thus our effective Rayleigh dissipation density is given by or, discretising as before, we have where we have once again absorbed constants associated with the discretisation of space into the parameters c, which are not to be confused with the parameters in the low-energy effective Hamiltonian. the phase relationships between degenerate modal coordinate Q and velocityQ, and the elastic τ elastic and damping τ damping forces. The diagram is drawn relative to Q, with a clockwise rotation indicating phase lag and counterclockwise indicating phase lead. B) steady-state entrainment to sinusoidal muscular forcing of magnitude A and frequency ω d leads to a constant neuromuscular-mechanical phase shift φ as determined from the mechanical transfer function. The force can always be decomposed into position (Q) and velocity (Q) components. The velocity component always completely cancels the damping force (otherwise oscillation would grow or decay) while the position component leads to an effective modification of the elastic force. C) the behaviour of the entrained system is equivalent to the behaviour of a rescaled, unforced system that experiences no damping and has shifted effective elastic properties.

B Incorporating neuromuscular forcing effects into the phenomenological parameters
We believe that the effective theory we have derived above from simple symmetry and naturalness arguments should describe not only the physics of the body, but should also capture the physics that arises due to the interaction of the neuromuscular system with the body.
Indeed, we believe that the general symmetry considerations we used to derive the Hamiltonian and Rayleigh function should be approximately true of the neuromuscular system. For instance, we do not believe that the action of the nervous system should be altered by the larva's overall position or orientation in space in the absence of interactions with the environment (including the absence of gravity, magnetic fields, temperature gradients, etc), and we do not believe the neuromuscular system should have an overall "handedness" given the bilateral symmetry of the musculature.
Furthermore, we very often observe in the construction of effective theories that one can ignore the complicated phenomena that occur on short length and time scales if we are willing to accept that our effective theory is only valid on long length and time scales. This is in fact one of the strengths of the effective theory approach -it forces us to be explicit about the domain of validity of the theory. For instance, disregarding the complexity of the DNA double helix, with all of its electrochemical interactions and beautiful detailed structure, and simply modelling the molecule as a uniform cylindrical rod, gives an excellent description of its physical behaviour on long length scales. This description breaks down when considering short-range effects, when aspects such as the chirality of the double helix become important [54]. Indeed, our results indicate that our effective theory of the larva correctly predicts the statistics and average trajectories during several larval behaviours, without having to make any real reference to the specifics of the nervous system, even though the real behaviour most certainly arises due to the interaction of the body with the neuromuscular system. This does strongly suggest that the effects of the neuromuscular system can be accounted for within our long-wavelength, low frequency theory.
In this section we attempt to provide a slightly more detailed sketch of why the phenomenological parameters in our theory, combined with some constraints related to the mechanical invariants, are able to account for neuromuscular effects. We will not aim to provide any strictly rigorous proof -the fact that our effective theory is so predictive is frankly surprising, and developing a precise understanding of why, for instance, equilibrium statistical mechanics is able to correctly predict the statistics of real larval behaviours, is an open problem.
We'll consider a single mode, described using the Hamiltonian and Rayleigh functions of our effective theory, then derive the mode's dynamics and subject it to some neuromuscular forcing. We'll apply a set of simplifying assumptions similar to those made in the construction of the effective theory (so that everything stays consistent), and show that we are left with forced dynamics that are equivalent to the original dynamics but with rescaled phenomenological parameters.
To begin, we write the Hamiltonian and Rayleigh function for a single mode which can be understood from Equations 7 and 8. The driven dynamics is given by the dissipative Hamiltons equationsQ we'll convert this to the familiar second-order form by subsituting the first equation into the second and rearranging to findQ where τ is the external force. In general, τ may be a function of the mechanical state variables Q,Q and time t. We'll make the dependence of the force on the mechanical state variables explicit (as though the force is being generated by a reflex, for instance) but the main ideas hold even if the relationship between the force and the mechanical state variables arises due to entrainment (for instance if the force is periodic, e.g. generated by a central pattern generator). We'll consider this second case later.
Next, we'll use Reynold's decomposition (which arises in the study of turbulent fluids) to split the force τ into the sum of an average component, that varies over the state space but does not depend on the time, and fluctuations around this average So far we have made no approximations, we have just written the force in a new, completely equivalent way. For our first approximation, we will assume that the correlation between the applied force and the mechanical state variables is entirely captured by the first, averaged, term, and that the fluctuations in the force are uncorrelated with the mechanical state variables. In this case, we can approximate the second term with a noise source ξ(t) - Next, we note that our effective theory is a free-field theory, with linear dynamics. It seems unfair to include a nonlinear external forcing term in our linear theory, given that we removed all possibility of nonlinear effective physics by restricting our attention to long length and time scales. Thus, we linearise the averaged force by expanding it in a Taylor series, where all of the bracketed expressions must be evaluated at Q =Q = 0, so they are simply constants. If we keep the first, constant term we will effectively shift the static equilibrium position of the mode. We use this kind of external force to usefully model self-righting behaviour in the main paper, but it is not what we are interested in now. Instead, let's see what happens when we include the next two, linear, terms, discarding the constant term and all higher terms. In this case we can write out the force τ as Substituting this expression into our dynamics gives which can be rearranged to giveQ Clearly c 2 and c 3 can simply be absorbed into the definitions of our phenomenological parameters ζ and ω, which gives us our original modal dynamics plus an external noise source -this is just a Langevin dynamics equation for the mode, which should give rise to equilibrium statistics under an appropriate choice of the noise ξ.
Thus, under simplifying assumptions of (1) decomposition of the applied neuromuscular force into a part that is correlated with the mechanical state and a part that is completely uncorrelated, (2) linearisation of the correlated part, and (3) an appropriate choice of noise term for the uncorrelated part, we arrive at a Langevin dynamics equation for the mode that should give rise to equilibrium statistics with shifted phenomenological parameters.
To show this, we started off by starting with a force τ that was explicitly dependent upon the coordinate Q and velocityQ as well as time. This may suggest that our line of reasoning is only valid for reflexive neuromuscular forces in which the coordinate and velocity are sensed and fed back, via the nervous system, to produce muscular forces. In fact, we believe that our reasoning should hold for any force that can be decomposed into a part that is correlated with the mechanical state plus an uncorrelated part. We'll demonstrate this by starting instead from a periodic external force τ that only explicitly depends upon the time and not the modal coordinate or velocity. In this case we can expand τ in a Fourier series where ω d is the fundamental frequency of the external driving force. The first term c 0 is just a constant, and as before it will function only to shift the modal equilibrium, so we discard it. Since we are interested in long length and time scales we keep only the c 1 and c 2 terms, i.e. we approximate the forcing by its fundamental component and throw away higher harmonics. This is in line with the modelling assumptions that went into building our effective theory, and is also supported by experimental evidence, since the fundamental component dominates the spectrum of the axial modes during peristalsis (if the forcing signal had important harmonics we would expect to see them in the modal spectra). In fact, if we choose the origin of our time axis appropriately we can keep just one of the remaining fundamental terms, so we have Our dynamics is thenQ The steady-state solution to this equation is well-known and is given by i.e. the modal coordinate oscillates at the forcing frequency, with a particular amplitude A and phase shift (relative to the forcing signal) φ. The amplitude A depends upon the forcing frequency ω d and the forcing amplitude c 1 while the phase φ depends only upon the forcing frequency ω d -the exact relationship is given by the transfer function H(ω d ) of the modal dynamics. The next few steps in our investigation of the periodic force can be made more precise by using the transfer function explicitly, however we believe this actually obscures the physical reasoning. Instead, we will use a closely related technique known as phasor analysis, in which we use vectors to pictorially represent the amplitude and phase relationships between quantities that are sinusoidally varying at a common frequency. This is done in Figure 10.
We start by plotting the coordinate Q as a vector of length A directed along the abscissa. We will measure all other quantities relative to this coordinate. The velocityQ always leads the coordinate by a 90 degree phase shift, so we plot it along the ordinate axis. The mode's elastic force is proportional to the modal coordinate Q, and directed opposite to it since τ elastic = −ω 2 Q, so this force is plotted in a negative direction on the abscissa, while the damping force is proportional to the modal velocityQ, and directed opposite to it since τ damping − 2ζωQ, so this force is plotted in a negative direction on the ordinate axis. This provides a complete picture of the unforced system ( Figure 10A).
Next we plot the external force τ at an angle −φ relative to Q ( Figure 10B). It should now be clear that we can decompose the force τ into a component that is in-phase with the velocityQ and a component that is in-phase with the coordinate Q. What are the magnitudes of these forces? Since we are in steady-state, we know that there can be no average power transfer through the mechanical system when we take into account both power flowing in to the system from the external force and power flowing out of the system due to the damping force. This tells us that the velocity-component of the external force must have exactly the same magnitude as the damping force -otherwise there would be an average transfer of power either into or out of the system; the oscillation would grow or shrink, and we would not be in steady-state. Thus, the effective damping force must be zero. Furthermore, we know that the frequency of the modal vibration must equal that of the forcing frequency. This tells us that the position-component of the force must sum with the elastic force in such a way as to give an effective elastic force −ω 2 d Q. These observations can be made precise by using the transfer function. We plot the force and its decomposition in Figure 10B. Ultimately, this is equivalent to a new unforced system that experiences no damping and has a shifted natural frequency ( Figure 10C), i.e. we have effectively made the transformation This can clearly be accounted for by a simple rescaling of the original effective parameters, ζ → 0 and ω → ω d . The fact that the amplitude of the output should stay constant in steady-state is readily accounted for by noting that our new system is derived from the effective Hamiltonian with no dissipation function, so this effective Hamiltonian (total energy) is conserved -this implies that the amplitude will neither grow nor decay.
To summarise, in this section we have provided two sketches of how the effects of neuromuscular forcing may be accounted for by the phenomenological parameters of our effective theroy. In the first, we assumed the external force could be written explicitly in terms of the coordinate, velocity, and time, but that it could be split into an averaged part that captures the correlation between the force and the mechanical state and a fluctuating part that is uncorrelated with the mechanical state; by linearising the correlated part and modelling the uncorrelated part as a noise source we arrived at a modified Langevin form of our original effective theory in which the phenomenological damping and elastic parameters were shifted and the mode was subjected to noisy forcing. In our second sketch, we assumed the external force was a periodic function of time that did not explicitly depend upon the mechanical state. Nevertheless, steady-state entrainment to the force lead to a similar modification of the phenomenological elastic and damping parameters. The final important effect, that arises in both approaches, is the appearance of a constant component of the force that acts to shift the equilibrium configuration of the effective theory. We use such a force in the main paper to account for self-righting behaviour.

C Dissipative axial self-interactions and the forward/backward frequency ratio
One of the main limitations of our effective theory is its inability to explain the asymmetric propagation speeds of forward and backward peristaltic waves using a single set of phenomenological parameters -to account for the asymmetry within the confines of our free-field theory we have to resort to using one set of parameters for forward peristalsis and another for backward peristalsis. In this section we provide a sketch of how to account for this asymmetry by incorporating self-interaction terms. We do not aim to guess the "correct" interaction term, rather we aim to demonstrate that it is possible to produce a qualitatively correct forward/backward asymmetry by incorporating interactions into our effective theory.
Since conservative Hamiltonian mechanics is time-reversal symmetric, we cannot account for the forward/backward asymmetry using a conservative interaction. To see this, suppose we build a Hamiltonian theory with a forward propagating wave solution. Hamiltonian theories are time-reversal symmetric, so the time-reversed wave solution must also be a valid motion generated by the theory. Clearly the speed at which the reversed wave propagates will be equal to the speed of the original forward wave, so their associated frequencies must match. We must therefore account for the asymmetric forward/backward frequencies using a dissipative interaction.
We will proceed by writing the dynamics for a pair of axial modes with degenerate frequency. We then modify the dynamics by adding a fairly general perturbation that should lead to a forward/backward asymmetry. We simplify the perturbation a little, leading to a cubic dissipative interaction force between the two modes. This is a self-interaction of the axial field with itself, since both modes belong to the axial field. Finally we prove that the resulting nonlinear theory has solutions in which the forward and backward propagating waves can have equal magnitude but different frequencies.
We begin from the dynamics for the first pair of axial modes X 1 , X 2 , writing their degenerate frequency without subscripts, as ω, for notational simplicity. To clarify our argument we initially neglect dissipation A simple way to make the propagation speeds for forward and backward motions asymmetric is to modify the degenerate frequency by some monotonic function of the generalised angular velocity for the pair of modes,φ = X 1Ẋ2 − X 2Ẋ1 , since a wave travelling in one direction hasφ > 0 while in the opposite direction it must haveφ < 0. We write this modification for the first modal coordinate as We can rearrange this equation to show that the new dynamics is just the original dynamics subjected to an interaction,Ẍ To simplify, we linearise the force with respect to the angular velocity, by assuming the force is analytic and expanding it in a Taylor series we ignore all quadratic and higher terms, and discard the leading constant term since it must simply act to shift the effective frequency parameter ω. The factor in brackets is a constant, since it is evaluated atφ = 0, so we simply write it as c for simplicity. The dynamics of the first mode is noẅ where we have used the definition of the angular velocity. The right-hand-side can clearly be interpreted as a cubic interaction force between the two modes, parameterised by a coupling constant c.
We now prove that this nonlinear model does indeed give rise to asymmetric forward/backward frequencies.
We take the ansatz where A is the amplitude of the travelling wave and ω s is its frequency, which can be either positive or negative depending on whether the travelling wave is forward-or backward-propagating. The ansatz can be differentiated to give the velocitiesẊ and differentiated once more to give the accelerations Substituting into the dynamics then gives The bracketed term can be simplified by using the trigonometric identity cos 2 a + sin 2 a = 1, and the entire equation can then be rearranged to bring out a common factor of A cos ω s t, leading to Since this equation must be true for all times t, the expression in brackets must be equal to zero. This gives us the following quadratic equation for the frequency of the solution ω s in terms of the coupling constant c, amplitude A, and free-field natural frequency ω, The solutions for ω s can then be found from the quadratic formula, There are two distinct frequencies, one negative and one positive, with unequal magnitudes, as desired. We plot the individual forward and backward frequencies and the forward/backward frequency ratio as functions of the coupling strength c in Figure 11, for the case A = ω = 1. In this instance the approximate 1 : 2 frequency ratio observed in the real larva can be obtained by setting c = − √ 2/2, as can be checked by the interested reader.

D Weak axial-transverse interactions and the crawling/exploration frequency ratio
We next consider the frequency relationship between the axial and transverse modes. Since the axial and transverse fields do not interact in our theory, it cannot tell us their relative frequencies. However, we can make good progress towards this aim by incorporating only a small nonlinear coupling into our theory.
Indeed, we have previously analysed the coupling between axial and mediolateral transverse degrees of freedom [30] using a simplified Hamiltonian for the axial and mediolateral transverse motions of the larva's head where q and φ are dimensionless variables characterising the axial stretch and mediolateral transverse bending angle, p q and p φ are the canonical momenta conjugate to these variables, λ is the axial-transverse frequency ratio, and is a global amplitude scale which multiplies the dimensionless quantities q, φ to give the real axial stretch and mediolateral bending angle of the head.
In the limit of small amplitudes the axial and transverse degrees of freedom are uncoupled, as in the effective theory we present this paper. We are interested in small coupling corrections to the uncoupled behaviour, so we expand the head Hamiltonian in a perturbation series in the amplitude parameter To keep our focus on the simplest model demonstrating axial-transverse coupling, we will keep only the zero-order small-amplitude term H 0 and the first order correction H 1 , finding and Since this coupled system is still fairly difficult to analyse, we choose to focus only on the transfer of energy from axial to transverse degrees of freedom. To do so, we will set the axial motion q to a prescribed function of time. We choose for this purpose q = cos(ωt), which could represent either an unperturbed axial modal vibration or a first-order approximation of the Fourier series expansion of a more complicated periodic function. Since q and p q are now prescribed functions of time, the terms in p 2 q and q 2 in the Hamiltonian will not effect the dynamics and can be discarded, leaving us with In order to simplify this Hamiltonian further, we next take a canonical transformation to new phase space coordinates given by Φ = p φ/λ , P = −λφ, i.e. we scale and interchange the coordinates and momenta, which allows us to group all parameters in one term and interpret the axial-transverse interaction as a sinusoidal modulation of frequency or, converting to second-order form by differentiating the first equation with respect to time and substituting into the second, we can write the dynamics in momentum space as which is in the form of the Matthieu equation [67]. It is well known that this equation exhibits the phonemenon of parametric resonance, in which the passively stable equilibrium at Φ = 0 becomes unstable for certain values of and ω, giving rise to solutions which grow with time.
For infinitesimal parametric perturbations (i.e. → 0), resonance occurs when the axial forcing frequency ω is exactly an even multiple of the natural frequency λ (i.e. ω = 2λ, ω = 4λ, etc.). For larger perturbations ( > 0), resonance occurs for a larger spread of frequencies centred on the even multiples of the natural frequency [55,67].
In the presence of friction, larger perturbations are required to produce resonance, and the magnitude of the required perturbation grows with the forcing frequency, so that only the low-order resonances are practically accessible [55,67]. The most readily excited resonance is therefore the 2 : 1 axial-transverse resonance, followed by the 4 : 1 resonance, and so on. Therefore, this line of argument suggests that if the larva favours the transfer of energy between axial and transverse degrees of freedom, the transverse frequency should be at a low-order, even subharmonic of the axial frequency. For instance, given an axial characteristic frequency of ≈ 1 Hz, the strongest parametric driving of the transverse degrees of freedom should occur when the transverse frequency is ≈ 1 2 , as observed in the real larva.

E Analytical eigendecomposition of D 2,c
We can find both Φ a and Λ a analytically by noting that D 2 is a circulant matrix. Indeed, the i'th eigenvector of an arbitrary circulant matrix is given by [68] Φ a,i = 1 where we have used Φ a,i to denote the i'th column of the eigenvector matrix Φ a , and z i = e 2πij N −1 is the i'th element of the (N − 1)'th roots of unity, with j = √ −1 the imaginary unit. Using Euler's complex exponential formula the k'th element of the i'th axial mode shape may be written The real and complex parts of each vector can be considered as independent mode shapes, so that the modes thus come in pairs with identical spatial frequency, For an arbitrary (N − 1) × (N − 1) circulant matrix with entries the eigenvalue corresponding to the i'th eigenvector is given by [68] In the case of D 2 we have c 0 = 2 and c 1 = c N −2 = −1, so that this reduces to However, the N − 1'th roots of unity satisfy z N −2 i =z i , where the bar indicates the complex conjugate. Therefore, the real part of z i can be found by using Euler's complex exponential formula, yielding By using the trigonometric identity 2 − 2cos (x) = 2sin x 2 we may now calculate the natural frequency of the i'th axial mode Let us now turn our attention to the eigendecomposition of D 2,f . The eigenvalue problem for this matrix can be stated as where the subscript t is intended to indicate that we are looking for the eigenvalue-eigenvector pairs related to twisting motions, k indexes over the segments of the model and i indexes over the eigenvector-eigenvalue pairs. For notational clarity we will drop the subscript t and index i until later in this subsection, and write the components of Φ t,k,i as v k . We will define the index k to run from 0 at the tail to n + 1 at the head.
Using this notation, the eigenvalue problem may be written component-wise as with the "boundary conditions" and v n − v n+1 = λv n+1 (120) We note that for the boundary conditions to be satisfied in general we must have which motivates us to introduce new variables w k = v k+1 − v k in which the boundary conditions can be written simply as w 0 = w n = 0. We note that so that the component-wise eigenvalue problem can be re-written but v k = w k−1 + v k−1 , so Noting that w k−1 − w k−2 = λv k−1 , we rewrite this as which can be rearranged to give Introducing 2α = 2 + λ and shifting the index k by 1 gives the reccurrence relation Considering α as an indeterminate, and assuming w 1 = 1 (which we can always satisfy by normalising the eigenvector, provided that w 1 = 0), we may write where U k is the k'th Chebyshev polynomial of the second kind. Our boundary condition at the head then tells us so that we can determine α, and thus the eigenvalues λ, by finding the roots of the (n − 1)'th Chebyshev polynomial of the second kind. These roots are well known, and tell us so that the i'th eigenvalue is determined by the equation rearranging gives or, by a trigonometric identity, This gives us n−1 non-zero eigenvalues. There must clearly also be a zero eigenvalue corresponding to uniform rotation of all segments, with corresponding eigenvector given by v k = c with c an arbitrary constant. We can incorporate the zero eigenvector into our expression above by writing so that the torsional natural frequencies may be written and the associated damping ratios are Each eigenvalue may now be substituted into the recurrence relations above in order to determine its corresponding eigenvector. Carrying out this procedure gives where the prefactor 2/N normalises the eigenvectors to unit length.

G Analytical eigendecomposition of D 4
To complete the eigendecomposition of our low-energy effective Hamiltonian, we now focus on the final eigenvalue problem where the subscript b is intended to indicate that we are looking for the eigenvalue-eigenvector pairs related to bending motions, k indexes over the segments of the model and i indexes over the eigenvector-eigenvalue pairs. As in the previous section, for notational clarity we will drop the subscript b and index i until later in this subsection, and write the components of Φ b,k,i as v k .
In the previous sections we exploited special properties of the second difference matrices to exactly compute their eigenvalues and eigenvectors. In particular, we exploited the circulant character of the second difference matrix with periodic boundary conditions to find its eigenvalues and eigenvectors using essentially a discrete Fourier transform, while we exploited the fact that the second difference matrix with free boundaries has a low bandwidth (non-zero elements only on the main diagonal and the two diagonals to either side, i.e. a bandwidth of 1) in order to write a tractable low-order recurrence relation for the eigendecomposition.
The fourth difference matrix we are now confronted with does not share these special properties -it is not circulant, and it has a higher bandwidth of 2 (i.e. five diagonals of the matrix have non-zero elements). This means that the discrete Fourier approach cannot be applied, and the recurrence relation will be of 5'th order and thus does not correspond to any Chebyshev polynomial.
Thus, we explore approximations to the eigendecomposition of D 4 . We begin by noting that our matrix can be decomposed into a product of rectangular forward and backward second difference matrices, i.e.
where we have noted that the backward and forward second difference matrices are simply related by taking the transpose, i.e. D 2,+ = D T 2,− . Since the nonzero eigenvalues of AA T and A T A are equal for any matrix A, this means that the nonzero eigenvalues of D 4 = D 2,− D 2,+ must be equal to those of the product D 2,+ D 2,− , i.e. the eigenvalues of which is the (N − 2) × (N − 2) fourth difference matrix with fixed boundary conditions, which we denote D 4,fixed . So far we have made no approximation, we have simply decomposed our original fourth difference matrix and proved that its nonzero eigenvalues must be identical to those for a related fourth difference matrix with different boundary conditions and a reduced dimensionality. Writing our problem in this way clarifies the mathematical issues we are facing -we need to find the eigenvalues of a real, symmetric, pentadiagonal Toeplitz matrix. There is currently no closed-form solution to this problem. However, we can use an asymptotics approach to find an exact solution in the limit N → ∞.
To proceed, we approximate the (N − 2) × (N − 2) matrix D 4,fixed by its circulant counterpart, which we denote D 4,c , with elements Clearly these two matrices are related by where the matrix C is zero everywhere apart from the six entries at the bottom left and top right corners (we choose the symbol C to stand for "corners"). Let us consider the relative error of this approximation as a function of the number of points N in our finite difference approximation. Clearly the element-wise approximation error is so that using a suitable matrix norm we can define a scalar measure of the absolute error or the relative error For concreteness let us use the L 2,1 matrix norm, which is simply the sum of the Euclidean norms of the columns of a matrix. Since the matrix C only changes through inclusion of more zero elements as N grows, C is a constant in this norm. Meanwhile, incrementing N by 1 introduces an additional column into D 4,fixed , always with the same 5 nonzero elements, so that the matrix norm should increase linearly as N → ∞. Therefore the relative error of our approximation should behave as ≈ 1/N as N increases, suggesting that the eigenvalues of the two fourth difference matrices should coincide aymptotically.
Next we note that the circulant fourth difference matrix can be factored into a product of circulant second difference matrices, Performing eigendecomposition thus gives so that we may identify the eigenvalue matrix of the circulant fourth difference matrix with the product of eigenvalue matrices of the circulant second difference matrix Λ 4,c = Λ 2,c Λ 2,c . Since the eigenvalue matrices are by definition diagonal, this means that the i'th eigenvalue of the circulant fourth difference matrix is the square of the i'th eigenvalue of the circulant second difference matrix. Using our analytical solution for the i'th nonzero eigenvalue of the second difference matrix thus gives the approximation which is our final expression for the eigenvalues of the fourth difference matrix with free boundary conditions.
We can find the eigenvectors of D 4,f analytically. Let us start from the diagonalisation condition we then perform a change of coordinates according to Φ = D 1 Φ 4 with D 1 the backward first difference matrix. We invert the change of coordinates to give Φ 4 = D −1 1 Φ, and insert into the diagonalisation equation to find Thus we seek vectors Φ to diagonalise the central matrix product D −1,T 1 D 4,f D −1 1 . Once found, we will be able to obtain expressions for the eigenvectors of the original matrix via Φ 4 = D −1 1 Φ.
We now note that our central matrix product gives us an augmented form of the (N − 1) × (N − 1) second difference matrix with free boundary conditions, i.e.
with the block diagonal matrix A given by which must have eigenvectors [0, Φ 2 ] T , where Φ 2 is an (N −1)-component eigenvector of the (N −1)×(N −1) second difference matrix with free boundaries. Thus, the eigenvectors of the fourth difference matrix with free boundaries are given by Meanwhile, the inverse of the backward first difference matrix is given by the N × N summing matrix which can be used in conjunction with our analytical expressions for Φ 2 from the previous section to find the analytical eigenvectors of D 4 . Carrying out this procedure gives the non-normalised eigenvectors H Detailed construction of the self-righting model We construct our model of self-righting by first finding the torsional equilibria of the larva. Since our effective theory is invariant under rotations, there exists a family of passive equilibria related by rotation about the axis of the larva. We can see that this is the case by starting from the effective Hamiltonian for torsional motion, or expressed in modal coordinates, where S i is the momentum conjugate to the torsional modal coordinate Θ i .
We then incorporate the effects of viscous power losses on the torsional motion using the effective Rayleigh dissipation function for torsional motion or in modal coordinates To find the passive mechanical equilibria we setθ = L =L = 0. Imposing this condition sets the Rayleigh dissipation function to zero, and eliminates the kinetic energy from the Hamiltonian The dissipative Hamilton's equation for the dynamics of the torsional angular momentum then giveṡ i.e. the torsional equilibria belong in the nullspace of D 2,f . Since this matrix is of rank N − 1, there must be a one-dimensional family of equilibria.
Inspection of the second difference matrix shows that these equilibria must correspond to θ = [c, c, · · · , c] T , i.e. at equilibrium, all segments are aligned with a constant rotation relative to the substrate, in line with our earlier symmetry considerations. Therefore, if we assume that one segment is fixed relative to the substrate (e.g. if the mouth hooks are attached to the substrate at the head) then we must have an equilibrium in which all body segments are aligned with the substrate. In other words, the attachment of one segment sets the particular value c in the equilibrium vector.
Since the Rayleigh dissipation function is negative definite in the absence of external forces, any passive motion must cause the body to lose energy and tend towards it's mechanical equilibrium at the minimum of the potential energy. Thus, regardless of starting configuration, if one segment is fixed relative to the substrate, then the rest of the body will eventually align with this segment and come to rest. For the larva, this means that attachment of the mouth hooks to the substrate is sufficient for self-righting.
We next ask how the body gets from its starting configuration to the righted, aligned configuration. We can find an answer by analysing the dynamics for the i'th torsional modë Θ i + ζ t,i ω t,iΘi + ω 2 t,i Θ i = 0 This second-order damped oscillator equations is well known, and the qualitative character of motion is entirely determined by the damping ratio ζ t,i . The i'th mode approaches the equilibrium configuration fastest in the case of critical damping for which ζ t,i = 1 and the approach to equilibrium occurs at a rate ζ t,i ω t,i = ζωλ t,i , while the equilibrium configuration is approached more slowly if ζ t,i is increased or decreased from this value.
Now we address the question of how the larva should move from some initial non-righted configuration into a configuration in which it may attach its mouth hooks to the substrate.
We assume that the absolute orientation of the head relative to the substrate θ h can be estimated by the larva, via some somatosensory channel. For instance, this may include a fusion of gravitometric, visual, and proprioceptive inputs. The current orientation of the head is clearly related to the superposition of the torsional modes. We argue that the greatest variance in the torsional motion should be along the first mode, which has the longest wavelength, lowest frequency, and lowest dissipation of all the torsional modes, and corresponds to a monotonic twist distributed across all segments of the animal.
Assuming the larva is prepared in a non-righted initial state in which all segments are at equilibrium at an orientation θ 0 , we may study how the larva may find a viable configuration for mouth hook attachment using a simple force input along the low-dissipation mode. To see this, we write the equation of motion for the modal coordinate asΘ + 2ζ t,1 ω t,1Θ + ω 2 t,1 Θ = Q There is clearly a linear relationship between Θ and θ h , with proportionality equal to the particular value of the first torsional mode vector at the head. We denote this constant of proportionality a, writing for the modal coordinate Θ = 1 a (θ h − θ 0 ). Substituting this expression into the equations of motion gives For mouth hook attachment, we expect that the application of the modal force Q should create a stable equilibrium with the head aligned with the substrate, so that θ h =θ h =θ h = 0. There are several choices of Q that could be used for this purpose -for instance, we could use PID control or a nonlinear controller. For simplicity we assume Q is a constant, as allowed for by our argument on the phenomenological effects of neuromuscular forcing (Section B), so that for the desired equilibrium we must have Indeed, substituting this expression back into the equations of motion gives θ h + 2ζ t,1 ω t,1θh + ω 2 t,1 θ h = 0 (166) so that the qualitative behaviour of θ h is again set by ζ t,1 , with the fastest possible self-righting occuring for critical damping when ζ t,1 = 1, for which the head will align with the substrate at a rate ζ t,1 ω t,1 exactly matching the rate at which the rest of the body segments will subsequently align with the head.

I Detailed construction of the peristalsis/rolling model
It is marked that the axial and transverse modes come in pairs with identical frequencies. Each such pair of modes is called a degenerate pair due to their relationship to the degenerate eigenvalues of the difference matrices in the effective Hamiltonian. Note that the axial modes in a degenerate pair have identical spatial frequencies, while those in the transverse case have completely identical spatial components.
The degeneracy of the eigenvalues of the difference matrices in the effective Hamiltonian suggests that the choice of corresponding eigenvectors is somewhat arbitrary. Indeed, in principal we may choose any linear combination of the two degenerate eigenvectors to form the basis for each degenerate pair. Even if we limit ourselves to considering only orthonormal degenerate bases, we are still free to rotate the degenerate basis vectors through any angle γ that leaves them in the same plane. Such a transformation cannot change the mechanics, and so the γ-rotations of the degenerate bases constitute additional continuous symmetries of the effective theory. A beautiful result in classical physics known as Noether's theorem tells us that these symmetries must be linked to additional conserved quantities. Let us make these ideas more concrete by considering a particular degenerate pair with modal coordinates X 1 , X 2 , conjugate momenta p 1 , p 2 , and natural frequency ω. To make the link between rotational symmetries and conservation laws most clear, we will briefly switch to using the Lagrangian framework of analytical mechanics. To do so, we begin with the Hamiltonian of the degenerate pair Taking the Legendre transform of the Hamiltonian yields the Lagrangian To make the rotational symmetry of this Lagrangian manifest, let us switch to polar coordinates A, γ defined via X 1 = A cos γ, X 2 = A sin γ. We will call A the degenerate amplitude and γ the degenerate phase. The Lagrangian is then so that the bracketed expression on the left hand side must be conserved. Indeed, we recognise this expression immediately as a generalised form of angular momentum, which we term the degenerate angular momentum. We now switch back to using the Hamiltonian framework by taking the Legendre transform of the Lagrangian, finding the degenerate Hamiltonian in polar coordinates where p = ∂L ∂Ȧ is the radial momentum conjugate to A and M is the degenerate angular momentum. This expression is simply a restatement of the total mechanical energy of the degenerate pair in polar coordinates. Note that we may now treat M as an arbitrary parameter and focus only on the dynamics of the degenerate amplitude A by introducing the effective potential energy U eff = M 2 A 2 + ω 2 A 2 . It should be clear that for a particular choice of degenerate angular momentum M there is now a minimum degenerate energy E = H(p, A) required for motion. Furthermore, the motion with this minimum energy corresponds to maintaining a constant degenerate amplitude, so that the trajectory of the system describes a circle in the original X 1 , X 2 configuration space. For finite energies exceeding this minimum, the motion is instead confined to an annulus in the X 1 , X 2 configuration space.
We now relate this picture back to the directly observable kinematics of the larva. First, we will convert from polar coordinates back to the degenerate modal coordinates X 1 , X 2 . To do so, we focus on the particular case in which E is set to its minimum value for a given M , corresponding to a conserved unit amplitude A = 1. We will set this value of M momentarily. Proceeding, we see that for our choice of A, M =γ is simply the angular velocity of the motion, so that γ = M t and the degenerate phase evolves linearly in time. Therefore, X 1 = cos M t and X 2 = sin M t. For this to coincide with our earlier results for the modal coordinates considered independently, we must have M =γ = ω, so that the modal coordinates X 1 and X 2 execute sinusoidal oscillations at the degenerate natural frequency ω with unit amplitude and a 90 • relative phase shift.
Next, we choose to interpret X 1 , X 2 as the two modal coordinates of the i'th axial degenerate pair with natural frequency ω a,i . Using our expressions for the i'th pair of axial eigenvectors, we can write the axial displacement x k of the k'th vertex of the midline as where we have dropped the normalising factor 1 √ N −1 . Using the identity cos (a) cos (b) ± sin (a) sin (b) = cos (a ± b), this further simplifies to Interpreting 0 ≤ k N −1 ≤ 1 as a spatial coordinate ranging over the undeformed configuration of the body, this is in the form of a sinusoidal travelling wave, and the choice of a minus or plus sign in the argument corresponds to the choice of a forward-or backward-propagating wave, respectively.
Restricting our attention to the axial degenerate pair with lowest frequency and lowest dissipation, and further assuming that the segments of the larva should be held in place without slipping when not moving, the translational speed of the larva should be Aγ = Aω a,1 .
Alternatively, we may choose to interpret X 1 , X 2 as the two modal coordinates of a transverse degenerate pair. For instance, the pair of transverse modes with the lowest frequency, and lowest dissipation, corresponds to C-bending in the mediolateral and dorsoventral planes. The rotational symmetry of the degenerate basis in this case corresponds to our arbitrary choice of mediolateral/dorsoventral axes, and the corresponding degenerate angular momentum conservation law gives rise to the rotational propagation of a C-bend around the body of the larva with angular velocity ω b,1 . Thus we see that there should be an exact correspondence between the angular velocity of rolling and the frequency of C-bending during unbiased behaviour, as we suggested in the main paper; this appears to be the case in the real animal.
In the presence of substrate constraints acting to hold the C-shaped bend in a fixed orientation in space (for The error function tends to 1 as its argument tends to ∞ and tends to −1 as its argument tends to −∞, so the definite integral becomes Substituting this into our expression for the partition function gives To calculate the remaining definite integral, we first calculate the antiderivative exp − β 2 Clearly our integral must diverge at A = 0 since the integrand includes a term proportional to A − 2, so to begin calculating our definite integral we instead find the limiting value of the antiderivative as A → 0 from above. We first distinguish two cases, corresponding to whether M is equal to zero or not. In the case M = 0 we have erf(a) → 0 and erf(b) → 0 as A → 0, so that the coefficients of the exponentials in the antiderivative are equal. However, the exponentials themselves are both equal to 1, since the exponents are 0, so that the exponentials cancel and the antiderivative must be equal to zero. Next considering the case M = 0 we have erf(a) → 1 and erf(b) → 1 as A → 0, so that the coefficients of the exponentials are both equal to zero. Thus we have lim Given that the value of the antiderivative vanished at our lower limit of integration, this last expression at our upper limit must be equal to the final definite integral in the partition function. Substituting this expression into the partition function then gives Substituting this into our expression for the partition function gives or, if we choose to measure phase space volumes in units of Planck's constant h, as is customary in classical statistical mechanics [49], we have where = h/2π is the reduced Planck's constant. This is our final expression for the partition function of a degenerate pair of modes subject to an average energy constraint and the constraint of the conserved momentum M = M 0 .