A physical theory of movement in small animals

16 All animal behaviour must ultimately be governed by physical laws. As a basis for understanding the physics 17 of behaviour in a simple system, we here develop an effective theory for the motion of the larval form of 18 the fruitfly Drosophila melanogaster, and compare it against a quantitative analysis of the real animal’s 19 behaviour. We first define a set of fields which quantify stretching, bending, and twisting along the larva’s 20 anteroposterior axis, and then perform a search in the space of possible theories that could govern the long21 wavelength physics of these fields, using a simplified approach inspired by the renormalisation group. Guided 22 by symmetry considerations and stability requirements, we arrive at a unique, analytically tractable free-field 23 theory with a minimum of free parameters. Unexpectedly, we are able to explain a wide-spectrum of features 24 of Drosophila larval behaviour by applying equilibrium statistical mechanics: our theory closely predicts the 25 animals’ postural modes (eigenmaggots), as well as distributions and trajectories in the postural mode 26 space across several behaviours, including peristaltic crawling, rolling, self-righting and unbiased substrate 27 exploration. We explain the low-dimensionality of postural dynamics via Boltzmann suppression of high 28 frequency modes, and also propose and experimentally test, novel predictions on the relationships between 29 different forms of body deformation and animal behaviour. We show that crawling and rolling are dominated 30 by similar symmetry properties, leading to identical dynamics/statistics in mode space, while rolling and 31 unbiased exploration have a common dominant timescale. Furthermore, we are able to demonstrate that 32 self-righting behaviour occurs continuously throughout substrate exploration, owing to the decoupling of 33 stretching, bending, and twisting at low energies. Together, our results demonstrate that relatively simple 34 effective physics can be used to explain and predict a wide range of animal behaviours. 35 1 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint

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, local Hamiltonian field theory (not illustrated in this figure).
larvae, we demonstrate that the statistics of several real larval behaviours can be surprisingly well described 84 using our equilibrium hypothesis. We explain the low-dimensionality of behaviour by Boltzmann suppression 85 of higher-energy, shorter-wavelength motions, with elementary low-energy, long-wavelength motions conse-86 quently dominating the animal's behaviour. We discuss the implications of these results for behavioural 87 control in Drosophila larvae and other animal systems.
Here Z is the partition function, which is defined so as to normalise the probability distribution. Because 128 we are dealing with continuous fields the partition function takes the form of a path integral, with the 129 measure Dψ intended to denote integration over all possible field configurations. We assume that the We further impose reflection symmetry; this is motivated by the bilateral symmetry of the larva, which we interpret as meaning that a left-handed deformation of the midline should have the same energy as a 145 right-handed deformation (so that the transformation y → −y should leave the Hamiltonian invariant).

146
Combining this left-right symmetry with overall rotational symmetry generates a greater range of reflection 147 symmetries -if we rotate the larva 180 • in the sagittal plane before performing a left-right reflection, we 148 have, in effect, performed a dorsal/ventral reflection, and we can clearly see that the Hamiltonian should where dependence on the derivatives of fields are implied, the term H I contains all quartic and higher order 155 terms, and the summation is over a set of purely quadratic, decoupled Hamiltonians, one for each field 156 x, y, z, θ ∈ ψ.

157
As we restrict attention to low energies (long wavelengths, low frequencies, and small field/derivative am-158 plitudes) both the higher order terms in H I and the higher derivatives within the quadratic Hamiltonians 159 H i become less relevant, and the Hamiltonian is dominated by the lowest order derivatives allowed by our 160 symmetry requirements (Appendix A). In this case we therefore arrive at Hamiltonian densities of fixed form H y = 1 2 (∂ t y) 2 + 1 2 c y (∂ 2 s y) 2 , H z = 1 2 (∂ t z) 2 + 1 2 c z (∂ 2 s z) 2 (4) where we have eliminated the coefficients of the first term in each Hamiltonian via a suitable rescaling of the 162 fields, in order to leave just four free parameters c i . We note in passing that these Hamiltonians are in the 163 same form as the Hamiltonian densities for stretching, bending, and torsion of a deformable rod in classical 164 linear elasticity theory, with stretching governed by the linear wave equation and bending governed by the 165 Euler-Bernoulli beam theory [43]. Our theory is thus strongly reminiscent of the phenomenological elasticity 166 models that have been used with great success in understanding mechanical properties of DNA and other 167 polymers [40,[44][45][46]; the key difference is our inclusion of kinetic energy, owing to the larva's relatively large 168 size and our interest in the larva's dynamics. 169 Having decoupled the Hamiltonian, the Boltzmann distribution and the partition function now factor into 170 separate contributions from each field with Z x = Dxe −βHx , and equivalent expressions for the other field-specific partition functions.
the Drosophila larva: T1 -T3 (thoracic), A1 -A8 (abdominal) 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 177 difference operators used in our discretisation. We impose vanishing force/moment (Neumann) boundary 178 conditions on the transverse and torsional fields, and periodic boundary conditions on the axial field. The

179
Neumann boundary conditions represent the absence of bending or twisting moments at the head and tail 180 extremities, while the periodic axial boundary conditions represent the presence of a visceral piston in the 181 larva which couples the motion of the head and tail. 182 We reduce the number of free parameters in our model by assuming that the transverse fields have an 183 "internal" circular rotational symmetry, i.e. we assume that the physics of mediolateral and dorsoventral 184 motions should be equivalent and write c y = c z , leaving three free parameters. This is motivated both by a 185 desire for simplicity in our theory, and by the approximately circular cross-section of the larva, which should 186 lead to a similar passive bending response in the mediolateral and dorsoventral directions. The introduction 187 of periodic boundary conditions to model the visceral piston of the larva introduces a similar "rotational" 188 symmetry into the effective theory. We discuss the implications of these symmetries later in the paper, 189 during our modal analysis and while studying peristalsis and rolling behaviours.

190
Our discretised Hamiltonian is then given by where we have once again combined the Hamiltonians for the separate fields, i.e. H = H x + H y + H z + H θ , 192 and we have indicated that the Hamiltonian corresponds to the total effective mechanical energy of the 193 midline, and is the sum of kinetic and potential energies.

194
In these expressions the lower-case bold symbols x, y, z, θ denote vectors containing measurements of the 195 scalar fields x, y, z, θ at uniformly sampled discrete points along the midline s i , while the vectors p x , p y , p z 196 contain the discrete measurements of the translational momenta and L contains discrete measurements of 197 angular momentum about the midline.

198
The parameters ω s , ω b , ω t measure the ratio of elastic to inertial effective forces and determine the abso-199 lute frequencies of stretching, bending, and twisting motions, respectively. Note these parameters are not Before exploring the behaviour of our effective theory, we first consider the effects of neuromuscular forcing 208 within the low-energy limit that we set above, since it is not intuitively clear whether our Hamiltonian 209 approach can effectively model these effects. To examine this issue we postulate that the Hamiltonian derived 210 above correctly encapsulates the low-energy dynamics, as well as the statistics, of our deformation fields. In 211 this case, we can derive partial differential equations governing evolution of the fields x, y, z, θ ∈ ψ by first 212 taking the Legendre transform of the Hamiltonian to find the Lagrangian density L(ψ) = i π i ∂ t ψ i − H(ψ) 213 and then finding the corresponding Euler-Lagrange equations which extremise the action S = dtdsL. Since 214 our Hamiltonian density is quadratic, so is the corresponding Lagrangian density, and the Euler-Lagrange 215 equation is therefore linear, and can be written with generality as 216 n,m (−1) n+m ∂ n s ∂ m t ∂L ∂(∂ n s ∂ m t φ) Here, φ ∈ ψ is used as a stand-in for any of the scalar fields x, y, z, θ, and P φ is a linear differential operator 217 encoding the Euler-Lagrange dynamics. In the presence of generalised forces acting on the field φ, the right 218 hand side is nonzero, and in general we will have instead represent an internal bending torque field, in which case we would have P µ = ∂ 2 s . Conversely, if φ represents 228 the axial field x(s, t) then the field τ (s, t) could represent internal tension forces, and we would have P µ = ∂ s .

229
The field τ may also have its own dynamics and may be reciprocally coupled to the mechanical field φ. 230 Assuming we remain within a linear, low-energy regime, we can write where P ν is constrained by symmetry. For instance, Galilean invariance will again set the lowest order of 232 derivative present within P ν , representing the idea that the neuromuscular system cannot produce forces 233 which depend upon the overall translation or rotation of the larva in space. However, in general P τ , repre-234 senting the dynamics of the τ field itself, is not required to satisfy any symmetries. 235 Now we note that we can apply the operator P µ to Equation 12, and assuming the commutativity of our 236 linear partial differential operators (i.e. assuming Clairaut's Theorem holds) we can write 237 P τ P µ τ = P µ P ν φ (13) but by Equation 11 we have P µ τ = P φ φ, so this is equivalent to This means the field τ can be eliminated by modifying the dynamics for the field φ. This can be rewritten 239 simply as Modal analysis of the effective theory predicts the principal components and dimensionality of larval behaviour. A) First six mode shapes predicted by our effective theory (coloured lines) compared to principal components of stretching (extracted from recordings of forward and backward peristaltic crawling behaviour), bending (extracted from recordings of unbiased substrate exploration, rolling, and self-righting behaviour) and twisting (extracted from recordings of self-righting behaviour) deformations. PCA was performed on individual larvae. Median principal component shown in black, interquartile range in dark grey, and 2nd-98th percentile range in light grey. Downward deflection corresponds to segmental compression (stretching modes), right-handed bending (bending modes) or clockwise rotation (twisting modes). Mode shape predictions are parameter-free. B) Predicted proportion of variance explained by each theoretical mode shape (coloured lines) compared to observed proportion of variance along each theoretical mode during real behaviour (median in black, interquartle range in dark grey, 2nd-98th percentile range in light grey). Data are shown for individual modes (top) and cumulatively (bottom). Larval behaviour is low-dimensional in that a large proportion of variance can be explained by the first few longest wavelength modes in each behaviour, as predicted by our theory. C) The frequency relationships observed between stretching modes during forward (i) and backward (ii) crawling, bending modes during unbiased exploration (iii), and twisting modes during self-righting (iii) compared to the fit by our model (boxplots show real data; coloured lines show model fit; fit obtained by tuning the three free parameters ω s , ω b , ω t using nonlinear least squares). (v) Frequencies of the first stretching mode during forward/backward crawling (fw/bk crawl) and the first "C-bending" mode during unbiased exploration (reorient) compared to the angular velocity of rolling. Our theory correctly predicts that the angular velocity of rolling should match the frequency of C-bending during exploration (p = 0.23, Mann-Whitney U-test) but incorrectly predicts that the frequency of forward and backward crawling should match (p = 0.009). 9 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Table 2. Number of modes required to explain 70-90% of variance during different larval behaviours, compared to model prediction. These percentages represent the range of commonly used thresholds when estimating the dimensionality of a system via PCA [49]. Reported values were computed by linear interpolation of data in Table 1. 10 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 mode shape describing the spatial pattern of movement, equivalent to a principal component vector -also 257 referred to as an eigenmaggot in the field [9] 258 The mode shapes are formally given by the eigenvectors of the difference matrices in our effective theory, so 259 that transforming to modal coordinates diagonalises these matrices. This means that the Hamiltonian can 260 be written in the form Here Q i denotes the i'th modal coordinate, P i denotes its conjugate momentum, and ω i is the characteristic 262 frequency of the mode. Frequency is written in terms of an absolute frequency scale ω, and the i'th eigenvalue 263 λ i of the relevant difference operator. The absolute frequency scale ω corresponds to the relevant free 264 parameter ω s , ω b , or ω t in our discretised Hamiltonian, depending upon which field the mode belongs to.

265
Clearly the ratio of frequencies for two modes of the same deformation field is independent of the absolute 266 frequency scale, depending only upon the eigenvalues of the difference matrices In terms of the modal coordinates, the Boltzmann distribution completely factors as where Z i = 2π/βω i is the partition function for the mode described by Q i , P i . Clearly, in this description, proportion of variance explained by each mode of a given deformation field is independent of the temperature 280 β and absolute frequency ω parameters, since we have This, therefore, enables us to make parameter-free predictions of both the mode shapes and the proportion 282 of variance explained by each mode (Table 1 and Figure 2; see Appendices E, F for calculations of the 283 eigenvalues and eigenvectors for axial and torsional deformations, and G for analytical approximations for

Principal components analysis
We compare the mode shapes predicted by our effective theory to those obtained via PCA of real larval 287 behaviour in Figure 2. There is a strikingly good agreement between theory and experiment across all of the 288 deformation modes and larval behaviours we studied : axial compression/expansion observed during forward 289 and backward peristaltic crawling, mediolateral transverse bending during unbiased substrate exploration 290 behaviour, mediolateral-dorsoventral transverse bending observed during self-righting and rolling behaviours, 291 and twisting observed during self-righting ( Figure 2A).

292
Furthermore, the proportion of variance explained by each of our predicted modes is surprisingly well ex-293 plained by equipartition of energy, with substantial variance explained by the first few long-wavelength modes 294 ( Figure 2B and have the same energy on average, but this amount of energy produces much larger amplitude excursions in 311 the long-wavelength modes than in their short-wavelength counterparts. As a result, observations of larval 312 behaviour are dominated by the long-wavelength modes.

313
One of the striking results of our analysis is that the eigenmaggot shapes are largely conserved between 314 unbiased behaviour and rolling behaviour -our theory also gives a strong explanation of why this must be 315 the case: these shapes are encoded in the effective physics of the midline. This, would otherwise be, highly 316 non-intuitive. For instance, it may be surprising that similar proportions of transverse variance are explained 317 by C-bending during unbiased behaviour (85.72%) and rolling (82.01%), despite this type of deformation 318 being more visually obvious during rolling. However, our equilibrium effective theory tells us that this should 319 be true simply because equipartition of energy amongst the modes should generically lead to a large amount 320 of variance being explained by C-bending.

321
The predictions so far have been independent of the particular values of our models' parameters. We can 322 now fix some of these parameters before moving forward, by considering the timescales of axial, transverse, 323 and torsional deformations during larval behaviour. In particular, we have tuned the three free natural 324 frequency parameters of our theory (which determine the relative timescales of the axial, transverse, and 325 torsional deformations, as well as the overall absolute timescale of our model) to fit the observed frequencies 326 for each mode. Indeed, we were able to find qualitatively good fits to these spectra ( Figure 2C), although 327 we were required to separately fit the spectrum of forward and backward peristaltic crawling as the overall 328 timescales of these behaviours differ from one another [21,50]. Explaining this discrepancy will require future 329 modification of our effective theory.
explain, this frequency ratio due to the lack of energetic interactions between axial and transverse modes.
However, in Appendix D we extend our theory to a higher energy regime by including a prototypical nonlinear 335 axial-transverse interaction (based on the geometric/kinematic nonlinearity in [19]); this interaction may lead 336 to transfer of energy between the axial and transverse fields via a 2 : 1 parametric resonance, thus predicting 337 an axial-transverse frequency ratio of ≈ 0.5, close to our measured value.

338
At this point, we are left with the temperature β as the only free parameter of the theory; this is joined by 339 an additional parameter in the following section. Although the Gaussian statistics of our theory successfully predict mode shapes and the proportion of 342 variance along each modal coordinate during real larval behaviour, they will in general fail to adequately 343 describe the probability distributions of these coordinates, which can be strongly non-Gaussian ( Figure 3;

344
Anderson-Darling normality test, p < 0.01 for all distributions of modal coordinates).

345
This discrepancy might be corrected by perturbatively modifying the Hamiltonian of our theory in order to 346 bring the predicted statistics closer to those we observe. Indeed, such perturbative modifications to a theory 347 are a cornerstone of the effective theory approach [35,36]. However, in our case we can make some progress 348 simply by looking closer at the symmetries of our Hamiltonian, without modification. In particular, we note words, the axial stretching/compression modes appear in pairs with identical frequency, and similarly for 355 every mediolateral bending mode there is a corresponding dorsoventral mode with the same frequency.

356
In our theory, due to Noether's theorem [51], the symmetries mentioned above are associated with mechanical 357 invariants: quantities which, like energy, do not change during the evolution of a closed system. The

358
derivation is presented below.

359
We start from the Hamiltonian for a pair of modes with degenerate frequency here, ω i is the common frequency of the two modes with coordinates Q 1 , Q 2 and conjugate momenta P 1 , 361 P 2 . To make the rotational symmetry of this Hamiltonian manifest, we take a canonical transformation to 362 polar coordinates in the Q 1 -Q 2 plane, by setting Q 1 = A cos φ, Q 2 = A sin φ, where A is the "amplitude", 363 or distance from the origin, in the plane, and φ is the "phase". In terms of these new coordinates we obtain 364 the Hamiltonian where p A is the canonical momentum associated with the amplitude coordinate A, and M is the canonical 366 momentum associated with the phase coordinate φ, corresponding to the angular momentum in the Q 1 -367 Q 2 plane. Note that φ does not appear in the Hamiltonian. This signifies that there is no energetic "cost" 368 13 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 associated with rotations in the Q 1 -Q 2 plane -this is our internal rotational symmetry. As a consequence, the angular momentum M is a mechanical invariant. This is a straightforward consequence of Hamilton's 370 equation for the dynamics of M , This can be interpreted as the usual Boltzmann statistics for H d multiplied by an additional Dirac delta 377 function which constrains M to take a fixed value M 0 . We calculate the partition function to be where the measure dΩ represents integration over phase space. Note that the most probable behaviour of the crawling the behaviour is similar, except the wave travels from head to tail and propagates at around half 392 the speed of its forward-propagating counterparts ( Figure 2). In both cases, the wave has a single "peak" of 393 maximal compression, but is spread across several segments, so that it appears largely as a broad sinusoid 394 of compression and expansion (recall that the first two axial modes alone explain around 55% of variance 395 during these behaviours, Table 1). In rolling, the larva enters a C-bend configuration (corresponding to 396 the first bending mode of our model -recall that this mode explains around 80% of variance during the 397 behaviour, Table 1). The C-bend stays in the plane of the substrate while the body rotates, so that the 398 C-bend propagates around the body axis in a clockwise or counterclockwise direction. For clockwise rolling, 399 the C-bend might first be directed to the right of the body, then ventrally, then to the left, then dorsally, 400 then back to the right; for counterclockwise rolling the sequence is reversed.

401
To the eye, rolling and peristalsis behaviours could scarcely look more different. However, owing to the 402 similar underlying rotational symmetry of the axial and transverse fields in our model, and the resulting 403 14 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; Fig 3. Similar physics underlies crawling and rolling behaviours (Previous page.) Ai) (top) A pair of modes with equal (degenerate) frequency ω can be described in terms of their respective modal coordinates (e.g. we use X 1 , X 2 for first pair of axial modes, and Y , Z for first mediolateral and first dorsoventral bending mode) or in terms of their amplitude A and phase θ. (bottom) Since our effective energy depends only on A and not θ, our model is symmetric under continuous rotations of the degenerate modes. Aii) This rotation corresponds to moving the peak of axial compression along the body, considering the head and tail as rigidly linked via the visceral piston [21] (rotation of X 1 , X 2 plane; top) or to rotating the direction of a C-bend around the body axis (rotation of Y , Z plane, bottom) Aiii) By Noether's theorem [47] the rotational invariance is linked to a conserved mechanical invariant (generalised momentum M ) corresponding to the momentum of a compression wave propagating along the body (top) or a C-bend rotating around the body axis (bottom). Aiv) Incorporating the conserved momentum into the theory as a constraint gives rise to an effective potential energy (U ef f , top) for the amplitude A [52]. The associated Boltzmann statistics for A can then be computed (middle); the most probable (average) trajectory minimizes the effective potential (orange lines) and has angular momentum M , corresponding to oscillation of the degenerate modes in quadrature (bottom) with amplitude M/ω and frequency ω. Av) The quadrature modal oscillations correspond to peristaltic compression waves (top; oscillation in X 1 , X 2 ) or rolling (bottom; oscillation in Y , Z). B, C) Fit of the model average trajectory and Boltzmann statistics to real data for forward crawling (B) and rolling (C); amplitudes were normalised so that M/ω = 1. i) Average trajectory (unit circle), real data (black), and a matched number of Markov chain Monte Carlo (MCMC) samples from the Boltzmann distribution (coloured circles) in the degenerate mode plane. ii) separate plots of the data in (i), including the analytical Boltzmann probability density function (pdf) iii) Representative trajectories of modal coordinates, amplitude, and phase over time (black) compared to average trajectory (coloured line) iv) Kernel density estimate (KDE, black) of the pdf of each quantity in (iii), computed over all trials and compared to predicted pdf (coloured line), with two predicted error bounds corresponding to either 2nd-98th percentile for MCMC samples drawn from the predicted pdf (light grey; see Methods) or interquartile range from the same MCMC samples (dark grey). v) cumulative distribution function (cdf) for each quantity in (iv) (same colorscheme; ), with results of Kolmogorov-Smirnov comparisons between the predicted cdf and the corresponding empirical distribution functions.
degeneracy of the axial and transverse spectra, both behaviours should share be described in mode space by 404 the degenerate statistics that we derived in the previous section, differing only in the values of the invariant 405 M 0 and the inverse temperature β.

406
In Figure 3 we plot experimental estimates for the first pair of degenerate modes during peristalsis (first 407 two axial compression modes, labelled X 1 , X 2 , Figure 3B) and during rolling (first mediolateral and first 408 dorsoventral bending mode, labelled Y , Z, Figure 3C), respectively, alongside a fit of our model to this 409 data. To perform our fit, we first eliminated M 0 by normalising the modes by the location of the peak of 410 the amplitude distribution (this peak lies at M 0 /ω in our model; Since the free parameter ω was fixed 411 in the previous section, normalisation is equivalent to fixing M 0 ). We then tuned β to fit the normalised 412 amplitude distribution using nonlinear least squares. We also plot our estimate of the (average) behavioural 413 trajectory both in the normalised mode space, where it corresponds to the unit circle, and as time series for 414 the individual modes, their amplitude A, and their phase φ. For both behaviours, the fit by our model is 415 qualitatively very good. Indeed, we find only one statistically significant difference between the predicted 416 and observed kinematic distributions, for the phase φ during rolling (KS test, p < 0.01, Figure 3Cv; all other 417 comparisons have p ≥ 0.11). This distribution is predicted to be uniform, but shows a peak near 0 • (which 418 corresponds to pure mediolateral C-bending), possibly reflecting the transitions to/from rolling at the initial 419 and final portions of our recordings. 420 We also note that our estimate of the average behavioural trajectory provides a novel experimental test of 421 our theory of rolling. In particular, our estimated trajectory tells us that the average angular velocity with 422 which the C-bend rotates about the body axis should simply be given by the frequency of the first transverse Kernel density estimate (black) of the probability density function (pdf) and cumulative density function (cdf) for the first five bending modes Y i (during unbiased exploration, left) and the first five twisting modes Θ i (during a mixture of unbiased exploration and self-righting, right) with Kolmogorov-Smirnov comparisons to the distributions predicted by our model. Each mode was normalised by the mode-specific standard deviation σ 2 i , which is fixed by our model up to an overall scaling by β. We determined this parameter by fitting the variance of the modes using nonlinear least squares. After normalisation, our model predicts each modal coordinate to have a zero-centred Gaussian distribution with unit variance (red and yellow lines). Although our model closely predicts the observed distributions (KS distance D ≤ 0.09 for all bending distributions and D ≤ 0.08 for 3 of 5 twisting distributions), all comparisons reached statistical significance (p < 0.01; see text for discussion).
all behaviours. Thus we predict that the average frequency of C-bending during unbiased behaviour should 425 match the average angular velocity of rolling. Indeed, we find no statistically significant difference between 426 our estimates of these two quantities ( Figure 2; two-sided Mann-Whitney U-test, p = 0.4). CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint with the configurational partition function Thus in the case M 0 = 0 the predicted modal distributions are again zero-centred Gaussians with variances 436 σ 2 i = 1/βω 2 b λ i , in line with our earlier results in which the momenta M were neglected.

437
We fit the free parameter β to fit the variances σ 2 i of the all the observed modal distributions during unbiased 438 behaviour. In Figure 4 we plot the first five modal distributions on axes normalised by our predicted

462
As for the bending deformations during unbiased exploration (see previous section), our Gaussian model 463 provides a qualitatively good approximation to the distributions of the torsional mode amplitudes (which we 464 were able to observe during a mixture of both self-righting and unbiased exploration behaviour; Figure 4).

465
However, all of the measured distributions show statistically significant differences when formally compared 466 to the Gaussian model predictions (KS test D = 0.16 ± 0.10, p < 0.01 for all comparisons). This is 467 unsurprising, given that we expect large torsional deformations to lie outside the domain of validity of our 468 low-energy theory. We believe our qualitative Gaussian approximation should still be seen as a useful starting 469 point for the development of higher energy theories, and should still provide a useful explanatory model.

470
18 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint Model-driven investigation of self-righting behaviour suggests the larva rights itself using twisting deformations A) In the twist-based model of self-righting (SR), the larva begins in an "upside-down" configuration (top), twists along its length to attach its mouth hooks to the substrate (middle) and then un-twists to align its body with the substrate (bottom). Alternatively, the larva may self-right via rolling (roll-based SR, see text). Bi) Representative trajectory showing the rotation angle of all body segments (apart from T2) about the body axis during the final portion of SR behaviour; note that the thoracic segments (T1, T3) remain in fixed alignment with the substrate throughout, while abdominal segments (A1-A8) rotate to align with them, until eventually the body is completely aligned with the substrate (yellow lines). Bii) overall rotation of the body (black) and estimated torsion along the body (grey, torsion estimated using first two torsional modes) for the trajectory shown in (Bi). C) overall rotation and torsion along the body are strongly negatively correlated (Pearson's r 2 = 0.81), with a slope near that predicted by the twist-based SR model (measured λ = −1.03, blue, vs predicted λ = −1, yellow), strongly suggesting Di, Dii) Kernel density estimates of the pdf and cdf for the first bending mode, measured in the plane of the substrate, during SR (black), substrate exploration (bright red), and rolling (dark red). The distribution during SR is closer to that for exploration (Kolmogorov-Smirnov distance D = 0.25) than that for rolling (Kolmogorov-Smirnov distance D = 0.32) and all distributions differ significantly (Kolmogorov-Smirnov test p < 0.01 for all comparisons) E) Furthermore, the angular velocity during rolling (dark red) matches the frequency of C-bending during unbiased exploration (bright red; p = 0.23, Mann-Whitney U test) while the angular velocity during self-righting differs from both (yellow, p < 0.01, Mann-Whitney U test), suggesting that self-righting is not driven by bending as in the roll-based SR model.

Model-driven experimental characterisation of larval self-righting behaviour 471
Finally, we turn our attention to the self-righting (SR) behaviour of the larva. This behaviour can be 472 observed by experimentally preparing larvae in an "upside-down" configuration with their dorsal surface 473 in contact with the substrate, after which the larva is able to turn itself so that its ventral surface is in 474 contact with the substrate, and the animals' normal spontaneous behaviour resumes [16]. At present, there 475 is no widely accepted theoretical, mechanistic explanation for how the larva accomplishes self-righting. An 476 obvious explanation would be that the larva simply uses the same mechanism that it uses to roll (as described 477 previously in this paper), which we will call the rolling-based SR model. However, an alternative explanation 478 has been posited [16,25,26], which suggests that during self-righting, and unlike during rolling behaviour, 479 19 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint the larva twists its body so as to attach its mouth hooks to the substrate, before rotating the rest of its segments to bring the body fully into alignment with the substrate ( Figure 5A). We will call this alternative explanation the twist-based SR model. Armed with our theory of 3-dimensional larval movement, we are 482 well placed to formalise these explanations and assess exactly how self-righting is accomplished. Since we 483 have already explored rolling behaviour earlier in this paper, we will now focus on the twist-based SR model. 484 We first observe that the passive mechanical equilibrium of our effective theory, i.e. the minimum energy 485 configuration, has all segments aligned with one another, since in this configuration there is no energy stored 486 in the twisting or bending of body segments relative to one another. Thus, after attachment of the mouth 487 hooks to the substrate, the minimum energy configuration (which corresponds to the most probable configu-488 ration due to the negative exponential weighting of energy in the Boltzmann distribution) will consist of all 489 segments being aligned with the substrate. In other words, attachment of the mouth hooks to the substrate 490 should be sufficient to produce full self-righting, via subsequent "passive" relaxation towards equilibrium. 491 We can formulate several experimental tests of this twist-based SR model. First, this model suggests that  Second, the overall (average) rotation of the larva relative to the substrate, captured by the zero-frequency 497 torsional mode, should be strongly correlated with the higher torsional modes. This follows because the 498 attachment of the mouth hooks to the substrate introduces a constraint on the torsional modes, which is 499 exactly what allows overall rotation of the larva to be driven by torsional deformations during self-righting.

500
In particular, the following identities must hold where θ head is the rotation of the head, v i,head is the element of the i'th torsional mode shape at the head, 502 and Θ i is the i'th torsional modal coordinate (with i = 0 corresponding to the zero-frequency overall rotation 503 mode). The overall rotation of the body at any instant is given by the average over segments θ i , which 504 corresponds exactly with the zero-frequency mode Θ 0 . This relationship can be interpreted as a simple 505 restatement of the fact that the overall rotation is proportional to the torsional deformation in the twist-506 based SR model. As our model explains, > 55% of the variance in torsional deformations is accounted for 507 by the first and second torsional modes alone. We should therefore have the approximation Experimentally, we have indeed observed a very strong correlation between the overall rotation and the (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; righted. This suggests that twist-based SR may be an ongoing process that occurs during other behaviours.
Thirdly, the rolling-based SR model suggests that self-righting should be driven by similar bending de-521 formations as those seen during rolling behaviour, rather than the twisting deformations predicted by the 522 twist-based SR model. Therefore, the rolling-based SR model predicts that the angular velocity associated with overall rotation during self-righting, i.e. d dt Θ 0 , should be equal to the average frequency of C-bending 524 during unbiased substrate exploration behaviour. This relationship holds during rolling, with the angular 525 velocity of rolling closely matching the average frequency of C-bending during exploration (median = 0.46 526 Hz for rolling; median = 0.43 Hz for exploration; Mann-Whitney U-test p = 0.23; Figure 5E, see Section 1.6).

527
However, during self-righting behaviour, the rate of overall rotation is lower than both the C-bending fre-528 quency during exploration, and the rate of overall rotation during rolling (median = 0.25 Hz for self-righting;

529
Mann-Whitney U-test p < 0.01 for both comparisons). This result is again in favour of the twist-based rather 530 than rolling-based model of self-righting.

531
Finally, we directly compare the distribution of C-bend amplitudes during self-righting, rolling, and unbiased 532 substrate exploration. We find that the C-bend amplitude distribution during self-righting is closer to that 533 observed during unbiased exploration than that observed during rolling (KS distance, D = 0.25 SR vs 534 unbiased exploration, D = 0.32 SR vs rolling; Figure 5D). However, all three distributions differ from one 535 another (KS test p < 0.01 for all comparisons): the C-bend distribution is unimodal and approximately 536 zero-centred Gaussian, the rolling distribution is bimodal (with peaks corresponding to large-amplitude left-537 handed or right-handed C-bending) and zero-centred, and the self-righting distribution is bimodal and not 538 zero-centred. Qualitatively, one of the peaks in the SR distribution appears to coincide with the peak of 539 the C-bend distribution for unbiased behaviour, and the other appears to coincide with one of the peaks of 540 the rolling distribution. Thus, although C-bending does not appear to play the same role during SR as it 541 does during rolling (since the first three experimental tests point to SR being driven by twisting rather than 542 bending), large-amplitude C-bending does appear to be a clear feature of SR. It is unclear why this should 543 be the case. One possibility is that C-bending during SR is a result of twist-bend coupling, a phenomenon 544 that occurs due to the geometrically nonlinear rod kinematics that become apparent at large deformations 545 (higher energy scales) [53], and which is not accounted for in our (linear) low-energy theory. Alternatively, 546 this observation may be due simply to our inability to measure segmental rotations greater than ∼ 90 • due 547 to occlusion of the anatomical landmarks used during tracking (see Methods), so that we are only able to 548 measure the final stages of SR; perhaps the initial stages of SR consist of rotation driven by C-bending, 549 as in rolling behaviour, and the larva only attaches its mouth hooks once it is close to being righted. This 550 would explain why we observe a large-amplitude peak in the C-bend distribution despite the fact that overall 551 rotation appears to be driven by twisting rather than bending. The resolution of this issue is unfortunately 552 beyond the scope of this paper; further progress must wait until either a higher energy theory is developed 553 or experimental limits are extended.

555
In this paper we develop an effective theory for the study of movement in small animals. In particular, our 556 work presents a model of the low-energy physics of the Drosophila larval AP axis (midline), and experimental 557 demonstration of this model's ability to predict real behavioural data. Our model is the first to describe the 558 fully 3-dimensional motions of the larval midline, including stretching along, bending of, and twisting about 559 the body axis. We believe it sets the stage for development of similar theories for the behaviour of a large 560 class of small animals with "slender" or "rod-like" body morphologies, including the invertebrate nematode 561 worm C. elegans and the vertebrate zebrafish D. rerio. Our theory may also be extended to the behaviour 562 of larger animals such as snakes.

563
The relation of theory to experimental observations 564 21 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint Our effective theory makes several predictions, many of which we have directly tested. For instance, our theory makes parameter-free predictions of the principal components of stretching, bending, and twisting, 566 as well as the proportion of variance explained by each component. These predictions closely match the 567 results of PCA applied to real larval deformations during a range of behaviours including forward/backward 568 peristaltic crawling, unbiased substrate exploration, rolling, and self-righting ( Figure 2). Since our effective 569 theory correctly predicts the proportion of variance associated with each mode, it may provide a simple 570 explanation for the low-dimensionality of larval behaviour, via Boltzmann suppression of short wavelength 571 modes. In this view, all of the deformation modes have the same average energy (due to equipartition), 572 but because the longer wavelength modes are "softer" than the "stiff" short wavelength modes, this energy 573 causes larger fluctuations in the longer wavelength modes. Thus, observations of behaviour are dominated 574 by the first few long wavelength modes.

575
Our theory is also able to predict low-dimensional average trajectories for the larva's rhythmic rolling and 576 peristalsis behaviours (Figure 3). It is also able to very closely predict the stretching and bending mode 577 distributions during these behaviours, and provides a good first approximation to the bending and twisting 578 mode distributions observed during substrate exploration and self-righting behaviours (Figure 3 and 4). As 579 a test of our model's utility, we were also able to use it to investigate the mechanism underlying self-righting 580 behaviour; our analysis in terms of our model's modes supported a view of self-righting being driven by 581 twisting deformations, rather than the bending deformations that drive rolling behaviour.

582
Given our model's explicit focus on the low-energy physics of the larva, we expect it to break down at the 583 higher energy scales governing short-wavelength and large amplitude deformations of the body. Indeed, 584 although our model provides a qualitatively very good fit to the mode shapes extracted from recordings of 585 real animal behaviour, the fit is certainly better for the dominant long-wavelength modes than for their short-586 wavelength counterparts (Figure 2). Furthermore, although we successfully predict the variance structure of 587 bending and twisting deformations, and our predicted Gaussian statistics for the corresponding distributions 588 in mode space are therefore qualitatively very good, our theory does not account for higher-order statistics. 589 We are thus unable to account for the heavy tails of these distributions, which represent relatively high-energy, 590 large-amplitude bending and twisting deformations. Given our theory's success at low energies, however, 591 we are hopeful that it may provide a useful starting point for a perturbative approach to understanding the 592 higher energy regime [35].

593
The model we propose raises several interesting conceptual questions -in particular, it remains unclear, 594 first, why our assumption of statistical equilibrium provides such a good description of larval behaviour; 595 and, secondly, why our theory works despite containing no detailed description of the larva's neuromuscular 596 system, environment, or detailed features of the (presumably nonlinear) biomechanics. We consider these 597 two questions in detail below.

598
Statistical equilibrium and its power to explain larval behaviour 599 We constructed our effective theory of the larval midline using an approach borrowed from statistical field  Crucially, our use of equilibrium statistical mechanics necessarily presupposes some mechanism that is not 612 22 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint contained directly in our effective theory. This should be clear since the effective theory is described by a quadratic Hamiltonian, which should give rise to linear deterministic dynamics [51], and so our theory 614 contains no mechanism by which trajectories may "wander" or "mix" in the phase space [54]; in other 615 words, the deterministic dynamics of our theory cannot explain its predicted statistics without modification.

616
This is actually a fundamental issue at the foundations of equilibrium statistical mechanics, and not specific 617 to our theory alone, since for very many systems to which equilibrium statistics are applied it is either not 618 possible to prove that the underlying dynamics of the system are capable of driving the system to equilibrium, 619 or it is possible to show that the dynamics are outright incapable of this feat (as in the case of our effective 620 theory) [40,54].

621
In some cases, effective equilibrium statistics may arise even though the underlying system is known to be out 622 of equilibrium [55,56]; presumably this must occur in the case of the larva, whose body mechanics experience 623 a strong flow of energy from the musculature. Indeed, it is interesting to consider how equilibrium statistics 624 may arise in larval behaviour. One argument is that small nonlinear modifications to our effective theory 625 should generically lead to energetic coupling between modes; it is quite possible that such coupling may lead 626 to chaotic dynamics, providing a mechanism by which the modes may thermalise. Our previous work on the 627 planar mechanics of unbiased exploration [19] suggests that beyond the low energy limit we have considered is present in larval behaviour, it may go some way towards explaining the animals' behavioural statistics.

634
Alternatively, the presence of noisy forcing may also be a factor in explaining thermalisation (we comment 635 on this below).

636
Accurate prediction of larval behaviour without a detailed description of the neuromuscular 637 system 638 We note that our theory is able to explain features of larval behaviour without including a detailed description 639 of the animals' nervous system. We argue that it is able to do so by capturing the essential combined effects (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint not mean that it has nothing to say about the role of the neuromuscular system. Indeed, by providing 663 a description in which the phenomenological parameters support multiple mechanistic interpretations (see 664 above), our theory is able to contribute to advancing our fundamental understanding of movement control, 665 by focusing effort on teasing apart the contributions due to intrinsic properties of body mechanics, from those 666 produced via neural actions and sensory processes. For example, the outcome of recent genetic screens aimed 667 at mapping genes (e.g. microRNAs) affecting larval movement [25,26,65] can now be pursued considering 668 effects of gene modulation and control on body mechanics, as well as, on neural development, function and 669 control. In addition, we can also envisage projections of our work into neural aspects of robotic design and 670 development: a theoretical understanding of what body mechanics can achieve per se, can help define the 671 minimal requirements needed to effectively command robotic movement.

672
More generally, our theory suggests that considerations of the effective physics of the body can go a long 673 way to explaining complex and diverse behaviours. In this context, the role of the nervous system may be 674 seen as being more modulatory in nature -shifting the global phenomenological parameters of the system -675 rather than representing a precise micro-management system acting to control the exact trajectories of the 676 animal body during behaviour.

677
Given the simplicity of our starting assumptions we believe our effective theory may be readily applicable   CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint      Larval midlines and contours were used to score bouts of rolling behaviour, as described elsewhere [71].

729
Larvae that were tracked for fewer than 5 seconds, or travelled less than one body length in distance, were (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint 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.
For calculation of larval deformations during self-righting, DeepLabCut (version 2.2b5) was used. Using a 738 standard DeepLabCut pipeline [73], at least 20 distinct frames were extracted per studied video using DLC's 739 K-means clustering. Up to 44 features were manually labelled per frame according to their visibility. These Posterior Spiracles). These manually annotated frames served as a training dataset fed into the ResNet-50 743 provided by DLC. The ResNet-50 was trained for 1,000,000 iterations via batch processing across NVIDIA 744 V100 GPUs. All video frames were then analysed using this trained network to produce estimated poses. We 745 then estimated the midline using the mean of the left and right contour markers, and estimated a tracheal 746 centroid by taking the mean of the left and right trachea markers. 747 We estimated the segmental rotation angles θ i between the substrate surface normal and the local sagittal 748 plane as θ i = tan −1 (d i /r i ), where d i is the measured distance between the midline and the tracheal centroid 749 at the i th segment, and r i is the radius of the body at that segment. We estimated the segmental diameter 750 (and hence radius) as being equal to the measured width of the larva in our camera projection, at each video 751 frame. in this dataset so did not need to be removed prior to PCA.

763
For comparison to our theory, the experimental data was also projected onto the modal basis of our effective 764 theory. We projected the abdominal segment length data onto the eigenvectors V 2,c of the circulant second 765 difference matrix D 2,c , which is identical to a spatial discrete Fourier transform basis (Appendix E), to 766 obtain estimates of the axial modal coordinates X i . We also normalised to correct for the truncation of 767 the axial modal basis (since only abdominal segments were measured). In particular, we modelled the axial 768 displacements as arising from a measurement process where P is the matrix which projects the "full" vector of axial displacements x real onto the subspace of 770 abdominal displacements (i.e. P models our lack of thoracic data). We wish to estimate the amplitudes X 771 of the axial modes. These mode amplitudes are defined via x = V 2,c X. Thus, we have 772 x measured = PV 2,c X real Ideally, we would like to invert this expression in order to find X real given x measured . We can use the 773 26 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint orthonormality of the eigenbasis to write V −1 2,c = V T 2,c (orthonormality follows from the fact that the circulant 774 second difference matrix is real and symmetric). However, the projection matrix is singular (as are all 775 projection matrices). We therefore use the Moore-Penrose pseudoinverse projection P T = P + to write an 776 initial estimate of the mode amplitudes using the truncated basis vectors V T 2,c P + , Using our expression for x real in terms of X real lets us write Which tells us our approximation will be improved by bringing the transformation on the left closer to the 779 identity transformation. We cannot rotate our basis vectors to orthogonalise them (since our purpose is 780 to estimate the deflections along the particular basis vectors in V 2,c ). However, we can at least bring our 781 transformation closer to an identity by normalising our new basis. To see this, imagine we have a pure 782 deflection X i along a single mode vector v i . Then our transformation would give the approximation which can clearly be made exact by normalising by the scalar v T i P + Pv i . Writing all such normalising factors 784 in the matrix N = diag(V T 2,c P + PV 2,c ) gives the improved approximation Finally, we again use our expression for X real in terms of x measured (along with the weak-inverse property of 786 the pseudoinverse P + PP + = P + ) to write our estimate of the mode amplitudes as In addition to applying this transformation, for figure 3 we further applied a brickwall (0.5-1.5Hz) bandpass 788 filter to the first pair of axial modes to remove noise and artifacts associated with the truncation.

789
The estimated abdominal and thoracic rotation angles θ were projected onto the eigenvectors of the free-790 boundary second difference matrix D 2,f . Since we were unable to measure the rotation of the second thoracic 791 segment T2, we estimated the torsional mode amplitudes using a similar method as for the axial modes, 792 applying the transformation where now V T 2,f is the eigenbasis of the free-boundary second difference matrix and P is the projection 794 matrix modelling our lack of measurements for T2.

795
The midline data from unbiased exploration, rolling, and self-righting behaviour experiments was projected 796 onto the eigenvectors of the free-boundary fourth difference matrix D 4,f . We estimated the eigenvectors 797 numerically using the numpy python numerics module [76] as the resulting eigenvectors are more accurate 798 27 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint than our analytical asymptotic approximations. 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    CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made   (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

33
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 In the main text, we summarised the construction of an effective Hamiltonian H for the axial, transverse, and torsional deformation fields x, y, z, θ of the larval midline. In this section we work through this construction 1013 in greater detail. Rather than starting directly from the Hamiltonian we will instead start from it's Legendre 1014 transform, the Lagrangian where ψ = [x, y, z, t] T is the vector of deformation fields and π is the conjugate canonical momentum vector.

1016
Assuming locality, we can write the Lagrangian as an integral L = dsL of a Lagrangian density. This As is common in field theory, we will occasionally refer to this Lagrangian density simply as the Lagrangian, 1020 despite the fact that this term technically refers to the integral of the density. Given a Lagrangian, the 1021 conservative dynamics of the deformation fields can be obtained from the Euler-Lagrange field equations 1022 n,m with ψ i ∈ ψ. Assuming the Lagrangian is analytic, we can perform a Taylor expansion and write Here we have adopted the convention that the zero'th order derivative is the identity, i.e. ∂ 0 ψ i = ψ i and 1024 summations are over all non-negative integer-order derivatives of the fields; indices on the constants b, c have 1025 been dropped (there is a unique constant for each summand).

1026
We can discard the constant a since it will not appear in the Euler-Lagrange equations and therefore does 1027 not effect the physics. Similarly, any linear terms (contained in the first summation) with g > 0 or h > 0 1028 will not contribute to the Euler-Lagrange equations and can also be discarded, leaving only linear terms of 1029 the form b i ψ i . Furthermore, only quadratic terms where the sum n + m + n + m is even will contribute to 1030 the Euler-Lagrange equations (as can be checked by the interested reader).

1031
We now apply our first symmetry requirement: that the Lagrangian be invariant under a reflection of each 1032 of the fields ψ i → −ψ i individually. This removes the remaining linear terms, and constrains ψ i = ψ j in the 1033 quadratic terms, i.e. to quadratic order the deformation fields are now completely decoupled. This lets us 1034 write the Lagrangian as a sum of field-specific quadratic Lagrangians plus an interaction part, where each of the first four L φ contains only the quadratic terms for the field φ ∈ {x, y, z, θ}, i.e.

34
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint while L I contains all higher order terms and interactions. We will now focus our attention on the quadratic Lagrangians, while continuing to use φ as a generic field symbol to denote any of the specific fields x, y, z, θ 1038 in our theory. Applying our second symmetry requirement: invariance under overall translations of the fields 1039 φ → φ + ∆φ. This eliminates the zero'th order derivative term where n = m = n = m = 0, corresponding 1040 to φ 2 , and we can also safely eliminate any other term containing a zero'th order derivative (since all such 1041 terms are completely equivalent to terms without a zero'th order derivative on going to the Euler-Lagrange 1042 equation). In other words, the only remaining terms must be of a minimum order with n + m ≥ 1 and 1043 n + m ≥ 1.

1044
Turning our attention to the Lagrangians for the transverse y and z fields, we can now apply our third 1045 symmetry requirement: invariance under overall rotations of the fields. For an infinitesimal rotation through 1046 angle δ about the center of the midline (s = 1/2), this takes the form φ → φ+(s−1/2)δ. This condition rules 1047 out the term c(∂ s φ) 2 , and thus all remaining terms must be of minimum order 2m + n ≥ 2 and 2m + n ≥ 2 1048 in L y and L z .

1049
Now we focus our attention on the low-energy (long-wavelength, low-frequency, low-amplitude) physics of 1050 the midline, by using a form of dimensional analysis inspired by the renormalisation group [36,40]. We 1051 first note that by focusing on low energies we are interested primarily in small deformations of the larva, 1052 characterised by small amplitudes of the fields and their derivatives, and we can therefore discard the 1053 interaction Lagrangian L I entirely (technically some terms in the interaction could be required to counteract 1054 instabilities of the fields [36], but we will soon see that we are lucky in this respect).

1055
We proceed by introducing a rescaling of space s = σs and time t = τ t into our theory, along with a 1056 compensatory rescaling of the field amplitude φ = φ . For τ > 1, σ > 1 we are essentially "zooming out" 1057 and viewing the larva on longer length and time scales. As we do so, the constants in our Lagrangians 1058 change, since on rescaling we have with the ratio of the constants before and after the rescaling transformation given by relative to the first two. What about the higher order terms? We can substitute = τ = σ into the growth equation (Equation 46) to find the relationship This tells us that in L x and L θ , all terms satisfying 2 < m + m + n + n , i.e. all terms above the lowest order, 1072 will shrink on rescaling. We can follow the same argument for the y and z fields. In this case we determined 1073 the minimum order terms to be given by the condition 2 = 2m + n and 2 = 2m + n . There are just two 1074 unique terms satisfying these conditions, namely the terms containing (∂ 2 s φ) 2 and (∂ t φ) 2 . Fixing these terms 1075 during rescaling gives = τ = σ 2 . The higher order terms in L y and L z then change on rescaling as which tells us that all terms satisfying 2 < n + n + 2m + 2m will shrink on rescaling; again, this corresponds 1077 to all terms containing derivatives higher than the minimal order.

1078
On repeated rescalings, most of the constants in our theory flow towards zero exponentially fast, and only 1079 the minimal order derivatives maintain non-zero coefficients. At the fixed point, we are therefore left with 1080 the Lagrangians where we have removed the coefficient of the first term in each Lagrangian through an appropriate redefinition for a single scalar field φ ∈ ψ = [x, y, z, θ] T .

1096
In what follows, we will perform a "momentum-space" coarse-graining. As a preliminary requirement, we 1097 must rewrite the field φ in terms of its spatial Fourier components φ k , which are labelled by their wavenumber 1098 k (which plays the role of momentum in quantum field theories, hence the nomenclature); intuitively, φ k 1099 measures the magnitude of a sinusoidal excitation of the φ field with wavenumber k. In terms of these 1100 Fourier modes, the partition function is given by where the high wavenumber (short wavelength) cut-off Λ is imposed because we expect our long-wavelength 1102 effective theory to break down beyond some short distance scale ∼ 1/Λ. To enact our coarse-graining, we 1103 introduce a new cut-off Λ < Λ, and redefine the Fourier modes φ k , the Hamiltonian H, and the partition 1104 function in terms of Λ . In particular, we split the modes φ k into low-wavenumber and high-wavenumber 1105 modes as where the low-wavenumber modes are given by and the high-wavenumber modes are given by A generic Hamiltonian can then be written as where H − , H + depend only upon either the low-or high-wavenumber modes, respectively, whereas H I 1110 governs interactions between the low-and high-wavenumber modes. The partition function now becomes The bracketed term on the right hand side can now be used to define a new effective, coarse-grained Hamil-In our particular case, we are interested only in the low-energy behaviour of the larval midline. We can therefore keep terms only up to quadratic order in H (see previous section) and neglect any higher-order 1116 contributions. In this case, the H must be diagonalised in Fourier space, so that H I = 0 and there are 1117 no interactions between the Fourier modes. This means that the integral within the brackets above must 1118 evaluate to a constant which is independent of the low-wavenumber modes, i.e. we can simply write so that the coarse-grained Hamiltonian is given by the expression Taking the logarithm and rearranging then gives In this section, we wish to provide an intuitive understanding for how a (fractional) quadratic Lagrangian 1132 density is able to account for the types of dissipation encountered within systems modelled by linear partial 1133 differential equations (PDEs). In particular, we will attempt to demonstrate how the general 2-dimensional 1134 scalar linear PDE of the form can be derived from a quadratic Lagrangian, via the Euler-Lagrange equation We start from the quadratic Lagrangian 1137 38 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 where q and φ are dimensionless variables characterising the axial stretch and mediolateral transverse bending 1156 angle, p q and p φ are the canonical momenta conjugate to these variables, λ is the axial-transverse frequency 1157 ratio, and is a global amplitude scale which multiplies the dimensionless quantities q, φ to give the real 1158 axial stretch and mediolateral bending angle of the head.

1159
In the limit of small amplitudes the axial and transverse degrees of freedom are uncoupled, as in the effective 1160 theory we present this paper. We are interested in small coupling corrections to the uncoupled behaviour, 1161 so we expand the head Hamiltonian in a perturbation series in the amplitude parameter 39 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 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 1167 from axial to transverse degrees of freedom. To do so, we will set the axial motion q to a prescribed function 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 1173 coordinates given by Φ = p φ/λ , P = −λφ, i.e. we scale and interchange the coordinates and momenta, which 1174 allows us to group all parameters in one term and interpret the axial-transverse interaction as a sinusoidal or, converting to second-order form by differentiating the first equation with respect to time and substituting 1179 into the second, we can write the dynamics in momentum space as 40 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 which is in the form of the Matthieu equation [78]. It is well known that this equation exhibits the phonemenon of parametric resonance, in which the passively stable equilibrium at Φ = 0 becomes unstable for 1182 certain values of and ω, giving rise to solutions which grow with time.

1187
In the presence of friction, larger perturbations are required to produce resonance, and the magnitude of the 1188 required perturbation grows with the forcing frequency, so that only the low-order resonances are practically 1189 accessible [52,78]. The most readily excited resonance is therefore the 2 : 1 axial-transverse resonance, We can find both Φ a and Λ a analytically by noting that D 2 is a circulant matrix. Indeed, the i'th eigenvector 1197 of an arbitrary circulant matrix is given by [79] 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 1199 i'th element of the (N − 1)'th roots of unity, with j = √ −1 the imaginary unit. Using Euler's complex 1200 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 1202 thus come in pairs with identical spatial frequency, For an arbitrary (N − 1) × (N − 1) circulant matrix with entries 41 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 the eigenvalue corresponding to the i'th eigenvector is given by [79] 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.

1207
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 can be stated as where the subscript t is intended to indicate that we are looking for the eigenvalue-eigenvector pairs related 1216 to twisting motions, k indexes over the segments of the model and i indexes over the eigenvector-eigenvalue 1217 pairs. For notational clarity we will drop the subscript t and index i until later in this subsection, and write v k+1 − 2v k + v k−1 = λv k , k = 1, · · · , n with the "boundary conditions" and 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 1224 written simply as w 0 = w n = 0. We note that but v k = w k−1 + v k−1 , so 1227 w k − w k−1 = λw k−1 + λv k−1 , k = 1, · · · , n 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 1230 w k+1 = 2αw k − w k−1 , k = 1, · · · , n − 1 Considering α as an indeterminate, and assuming w 1 = 1 (which we can always satisfy by normalising the 1231 eigenvector, provided that w 1 = 0), we may write 43 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made or, by a trigonometric identity, 1239 λ i = −4 sin 2 iπ 2n , i = 1, · · · , n − 1 This gives us n−1 non-zero eigenvalues. There must clearly also be a zero eigenvalue corresponding to uniform 1240 rotation of all segments, with corresponding eigenvector given by v k = c with c an arbitrary constant. We 1241 can incorporate the zero eigenvector into our expression above by writing 1242 λ i = −4 sin 2 (i − 1)π 2n , i = 1, · · · , n 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 corre-1245 sponding eigenvector. Carrying out this procedure gives 1246 44 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made where the prefactor 2/N normalises the eigenvectors to unit length.

G Analytical eigendecomposition of D 4 1248
To complete the eigendecomposition of our low-energy effective Hamiltonian, we now focus on the final 1249 eigenvalue problem where the subscript b is intended to indicate that we are looking for the eigenvalue-eigenvector pairs related 1251 to bending motions, k indexes over the segments of the model and i indexes over the eigenvector-eigenvalue 1252 pairs. As in the previous section, for notational clarity we will drop the subscript b and index i until later 1253 in this subsection, and write the components of Φ b,k,i as v k .

1254
In the previous sections we exploited special properties of the second difference matrices to exactly compute 1255 their eigenvalues and eigenvectors. In particular, we exploited the circulant character of the second difference 1256 matrix with periodic boundary conditions to find its eigenvalues and eigenvectors using essentially a discrete 1257 Fourier transform, while we exploited the fact that the second difference matrix with free boundaries has 1258 a low bandwidth (non-zero elements only on the main diagonal and the two diagonals to either side, i.e. a 1259 bandwidth of 1) in order to write a tractable low-order recurrence relation for the eigendecomposition.

1260
The fourth difference matrix we are now confronted with does not share these special properties -it is not 1261 circulant, and it has a higher bandwidth of 2 (i.e. five diagonals of the matrix have non-zero elements). This 1262 means that the discrete Fourier approach cannot be applied, and the recurrence relation will be of 5'th order 1263 and thus does not correspond to any Chebyshev polynomial.

1264
Thus, we explore approximations to the eigendecomposition of D 4 . We begin by noting that our matrix can 1265 be decomposed into a product of rectangular forward and backward second difference matrices, i.e.

1266
D 4 = D 2,− D 2,+ where we have noted that the backward and forward second difference matrices are simply related by taking 1267 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, 1268 45 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; i.e. the eigenvalues of use an asymptotics approach to find an exact solution in the limit N → ∞.

1277
To proceed, we approximate the (N − 2) × (N − 2) matrix D 4,fixed by its circulant counterpart, which we where the matrix C is zero everywhere apart from the six entries at the bottom left and top right corners 1281 (we choose the symbol C to stand for "corners"). Let us consider the relative error of this approximation 46 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 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 1287 columns of a matrix. Since the matrix C only changes through inclusion of more zero elements as N 1288 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 1290 as N → ∞. Therefore the relative error of our approximation should behave as ≈ 1/N as N increases,

1291
suggesting that the eigenvalues of the two fourth difference matrices should coincide aymptotically.

1292
Next we note that the circulant fourth difference matrix can be factored into a product of circulant second 1293 difference matrices, Performing eigendecomposition thus gives so that we may identify the eigenvalue matrix of the circulant fourth difference matrix with the product of which is our final expression for the eigenvalues of the fourth difference matrix with free boundary conditions.

1301
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. 1303 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 1305 able to obtain expressions for the eigenvectors of the original matrix via Φ 4 = D −1 1 Φ. 1306 We now note that our central matrix product gives us an augmented form of the (N − 1) × (N − 1) second 1307 difference matrix with free boundary conditions, i.e. 1308 47 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021.
with the block diagonal matrix A given by (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ; https://doi.org/10.1101/2020.08.25.266163 doi: bioRxiv preprint 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 H = 1 2 p 2 1 + p 2 2 + ω 2 X 2 1 + X 2 2 (124) Taking the Legendre transform of the Hamiltonian yields the Lagrangian 1333 L = 1 2 Ẋ 2 1 +Ẋ 2 2 − ω 2 X 2 1 + X 2 2 (125) To make the rotational symmetry of this Lagrangian manifest, let us switch to polar coordinates A, γ defined 1334 via X 1 = A cos γ, X 2 = A sin γ. We will call A the degenerate amplitude and γ the degenerate phase. The

1335
Lagrangian is then where p = ∂L ∂Ȧ is the radial momentum conjugate to A and M is the degenerate angular momentum. This 1343 expression is simply a restatement of the total mechanical energy of the degenerate pair in polar coordinates.

1344
Note that we may now treat M as an arbitrary parameter and focus only on the dynamics of the degenerate 1345 amplitude A by introducing the effective potential energy U eff = M 2 A 2 + ω 2 A 2 . It should be clear that 1346 for a particular choice of degenerate angular momentum M there is now a minimum degenerate energy 1347 E = H(p, A) required for motion. Furthermore, the motion with this minimum energy corresponds to 1348 maintaining a constant degenerate amplitude, so that the trajectory of the system describes a circle in the 1349 original X 1 , X 2 configuration space. For finite energies exceeding this minimum, the motion is instead 1350 confined to an annulus in the X 1 , X 2 configuration space. 1351 We now relate this picture back to the directly observable kinematics of the larva. First, we will convert from 1352 polar coordinates back to the degenerate modal coordinates X 1 , X 2 . To do so, we focus on the particular 1353 49 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020 case in which E is set to its minimum value for a given M , corresponding to a conserved unit amplitude is simply the angular velocity of the motion, so that γ = M t and the degenerate phase evolves linearly in 1356 time. Therefore, X 1 = cos M t and X 2 = sin M t. For this to coincide with our earlier results for the modal 1357 coordinates considered independently, we must have M =γ = ω, so that the modal coordinates X 1 and X 2 1358 execute sinusoidal oscillations at the degenerate natural frequency ω with unit amplitude and a 90 • relative 1359 phase shift.

1360
Next, we choose to interpret X 1 , X 2 as the two modal coordinates of the i'th axial degenerate pair with 1361 natural frequency ω a,i . Using our expressions for the i'th pair of axial eigenvectors, we can write the axial 1362 displacement x k of the k'th vertex of the midline as 1363 x k = cos (ω a,i t) cos 2πi k N − 1 ± sin (ω a,i t) sin 2πi k N − 1 where we have dropped the normalising factor 1 √ N −1 . Using the identity cos (a) cos (b) ± sin (a) sin (b) = 1364 cos (a ± b), this further simplifies to 1365 x k = cos ω a,i t ± 2πi k N − 1 Interpreting 0 ≤ k N −1 ≤ 1 as a spatial coordinate ranging over the undeformed configuration of the body, 1366 this is in the form of a sinusoidal travelling wave, and the choice of a minus or plus sign in the argument 1367 corresponds to the choice of a forward-or backward-propagating wave, respectively.

1368
Restricting our attention to the axial degenerate pair with lowest frequency and lowest dissipation, and 1369 further assuming that the segments of the larva should be held in place without slipping when not moving, 1370 the translational speed of the larva should be Aγ = Aω a,1 .

1371
Alternatively, we may choose to interpret X 1 , X 2 as the two modal coordinates of a transverse degenerate 1372 pair. For instance, the pair of transverse modes with the lowest frequency, and lowest dissipation, corresponds 1373 to C-bending in the mediolateral and dorsoventral planes. The rotational symmetry of the degenerate basis 1374 in this case corresponds to our arbitrary choice of mediolateral/dorsoventral axes, and the corresponding 1375 degenerate angular momentum conservation law gives rise to the rotational propagation of a C-bend around 1376 the body of the larva with angular velocity ω b,1 . Thus we see that there should be an exact correspondence 1377 between the angular velocity of rolling and the frequency of C-bending during unbiased behaviour, as we 1378 suggested in the main paper; this appears to be the case in the real animal.

1379
In the presence of substrate constraints acting to hold the C-shaped bend in a fixed orientation in space (for 1380 instance the body may be held parallel to the plane of the substrate due to either "repulsive" ground-reaction exp − β 2 M 2 A 2 + ω 2 A 2 dA = 1 2 π 2βω 2 e β|M |ω (erf(a) − 1) − e −β|M |ω (erf(b) − 1) with 1445 a = βM 2 + A 2 βω 2 √ 2A (154) and Clearly our integral must diverge at A = 0 since the integrand includes a term proportional to A − 2, so to 1447 begin calculating our definite integral we instead find the limiting value of the antiderivative as A → 0 from 1448 above. We first distinguish two cases, corresponding to whether M is equal to zero or not. In the case M = 0 Given that the value of the antiderivative vanished at our lower limit of integration, this last expression 54 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Next, we compute the Helmholtz free energy from the partition function by using Which gives 1467 F = log 2β 2 ω + β|M |ω β = log 2β 2 ω β + |M |ω (161) 55 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 16, 2021. ;https://doi.org/10.1101https://doi.org/10. /2020