Mechanics

To explore larval crawling and turning behaviours, we choose to describe the motion of the larval body axis (midline) in a plane parallel to the substrate (Fig 1, S1 Fig). The larval body is capable of more diverse motions including lifting/rearing [21], rolling [45], digging [46], self-righting / balancing, and denticle folding which we have recently observed to occur during peristalsis (S1 Video). However, while exploring flat surfaces, the larva displays fairly little out-of-plane motion (neither translation perpendicular to the substrate nor torsion around the body axis) and only small radial deformations [47]. Furthermore, the majority of ethological characterisations of larval exploration treat the animal as if it were executing purely planar motion [4, 6, 8–13, 48]. A planar model is thus a reasonable abstraction for the issues addressed in this paper, i.e., the generation of peristalsis, bending, and substrate exploration.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Our model of axial and transverse motion over a planar substrate.

The midline of the larva is modelled as a set of discrete point masses interacting with each other via linear, damped translational and torsional springs, and with the environment via Coulomb sliding friction. We model the larva’s incompressible coelomic fluid by constraining the total length of the midline to remain constant (see main text). Quantities used to describe deformations of the body, and interaction with the substrate, are shown in S1 Fig.

https://doi.org/10.1371/journal.pcbi.1006635.g001

The segmented anatomy of the Drosophila larva allows us to focus our description of the midline to a set of N = 12 points in the cuticle, located at the boundaries between body segments and at the head and tail extremities. We assign each point an identical mass, and measure its position and velocity relative to a two dimensional cartesian coordinate frame fixed in the substrate (the laboratory or lab frame). We therefore have N_DOF = 2N = 24 mechanical degrees of freedom. We note that our assumption of a uniform mass distribution along the midline is somewhat inaccurate, since thoracic segments are smaller than abdominal segments. However, simulations with non-uniform mass distribution give results which are qualitatively close to those presented here.

We assume that the larval body stores elastic energy in both axial compression/expansion and transverse bending, due to the presence of elastic proteins in the soft cuticle. We assume that energy is lost during motion due to viscous friction within the larva’s tissues and sliding friction between the body and the substrate. Sliding friction also allows shape changes (deformations) of the body to cause motion of the larva as a whole relative to the substrate (centre of mass motion).

Since the mechanical response of the larva’s tissues is yet to be experimentally determined, we assume a linear viscoelastic model. This is equivalent to placing linear (Hookean) translational and torsional springs in parallel with linear (Newtonian) dampers between the masses in the model, as shown in Fig 1, or to taking quadratic approximations to the elastic potential energy and viscous power loss (as in S1 Appendix). We note that the accuracy of the approximation may decrease for large deformations, in which nonlinear viscoelastic effects may become important.

As with larval tissue mechanics, there has been little experimental investigation of the forces acting between the larva and its environment. We therefore assume a simple anisotropic Coulomb sliding friction model, in which the magnitude of friction is independent of the speed of motion, but may in principle depend upon the direction of travel. This anisotropy could be thought of as representing the biased alignment of the larva’s denticle bands, or directional differences in vertical lifting or denticle folding motions which are not captured by our planar model. A mathematical formulation of our sliding friction model is given in S1 Appendix.

In addition to power losses due to friction, we also allow power flow due to muscle activation. For the sake of simplicity, we choose to allow only laterally symmetric muscle tensions. In this case, the musculature cannot directly cause bending of the midline, and can only explicitly drive axial motions. We will see later that even indirect driving of bending motion can lead to surprisingly complex behaviour, due to energetic coupling of axial and transverse degrees of freedom.

The choice to neglect asymmetric muscle tensions is made in order to simplify our model and provide a clearer illustration of the potential role of body mechanics in generating complex larval behaviour. We note that there is only one way for muscle activations to be symmetric—if we were to allow asymmetry we would have to specify exactly what form that asymmetry should take, and we have little empirical or theoretical grounds on which to do so. Nevertheless, there are some interesting cases which may be considered in passing—the presence of a constant torque about the model’s segment boundaries should cause a shift in the equilibrium posture towards a resting curved shape. The presence of torques which are a linear function of the local body bending angle or local angular velocity should shift the effective transverse stiffness or viscosity of the body S5 Appendix. In this sense the model presented here could be considered to already include the effect of asymmetric muscle tensions, they have simply been incorporated into the passive stiffness and viscosity of the body. We have recently developed an extension of the model presented here which uses a similar local reflex to modulate the body’s effective transverse viscosity in proportion to a stimulus input, allowing the model to exhibit taxis behaviour [49].

Finally, we model the internal coelomic fluid of the larva. Given the extremely small speed of the fluid motion compared to any reasonable approximation to the speed of sound in larval coelomic fluid, we can safely approximate the fluid flow as incompressible [50]. This would ordinarily require that the volume contained within the larval cuticle remain constant. However, since we are modelling only the motion of the midline and neglecting radial deformations, we constrain the total length of the larva to remain constant. We note that this constraint is not entirely accurate to the larva, as the total length of the animal has been observed to vary during locomotion [47]. Nevertheless, for the sake of simplicity we will continue with this constraint in place, noting that this approximation has been used with success in previous work focused on peristalsis [40, 43], and that there is experimental support for kinematic coupling via the internal fluid of the larva [3]. We note that we satisfy the incompressibility condition only approximately in some sections (Model behaviour—Conservative chaos, Dissipative chaotic deformations, and Deterministic exploration), by introducing an additional potential energy associated with the constraint, which produces an energetic barrier preventing large changes in the total length of the body (see S1 Appendix for details of this approximation along with specifics of the mathematical formulation of our mechanical model).

Note that in the absence of transverse bending, the total length constraint causes the head and tail extremities of the larva to become mechanically coupled and move in unison [40, 43]. The axial mechanics thus has periodic boundary conditions, and the most anterior (T1) and posterior (A8) segments of the larva may be considered adjacent. This means, for instance, that a compression wave travelling from tail to head will cause motion of the tail on termination at the head, thus initiating a new compression wave. This view also allows us to reason about what should happen if we relax the total length constraint. In particular, if we were to replace the direct coupling of head and tail by a viscoelastic coupling, representing the capacity for storage and dissipation of energy within the internal fluid or in radial expansion of the cuticle, the axial mechanics would still have periodic boundary conditions but would now have a step change in mechanical impedance. Waves hitting such “sudden” impedance boundaries in their transmission media will generally be partially transmitted (i.e. passing directly from head to tail in the larva) and partially reflected (i.e. changing direction and moving backwards from the head extremity), providing one possible cause of transitions between forward and backward locomotion in the animal. As will be seen, however, these transitions may occur even in the absence of an impedance discontinuity, and we will continue with the total length constraint in place in order to simplify our model.

Neuromuscular system

Let us now consider how we should use muscle activity to produce locomotion. There are two basic requirements. First, we must have power flow into the body from the musculature, so that the effects of friction may be overcome and the larva will not tend towards its equilibrium configuration. Second, we must be able to produce a net force on the centre of mass of the larva, so that it can accelerate as a whole relative to the lab frame. Note that in this section, we motivate the neural circuits in the model from this purely functional point of view, but will present relevant biological evidence in the discussion.

To satisfy the first criterion, let us examine the flow of power into the body due to the action of the musculature (1)

Here, q_i describes the change in length of the i’th body segment away from its equilibrium length, is the rate of expansion of the i’th body segment, b_i is a (positive) gain parameter, MF_i is a (positive) dimensionless control variable representing muscle activation, and the product b_iMF_i is the total axial tension across the i’th body segment. From this expression, it is clear that if we produce muscle tensions (MF_i > 0) only while segments are shortening (), we will always have positive power flow into the body (P > 0). This is a mathematical statement of the requirement for the larva’s muscles to function as motors during locomotion, rather than as springs, brakes, or struts [33].

A simple way to fulfil this condition is to introduce a segmentally localised reflex circuit (Fig 2, [40]). We place a single sensory neuron in each segment which activates when that segment is compressing (). Each sensory neuron then projects an excitatory connection onto a local motor neuron, which in turn projects to a muscle fibre within the same segment. Assuming for now that there are no other influences on the motor neurons, so that sensory activation implies local motor neuron activation, segmental shortening will produce an immediate muscle tension serving to amplify compression of the segment and thus counteract frictive energy losses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The neuromuscular model.

A local reflex amplifies motion via positive feedback: sensory neurons SN activate during segmental shortening, exciting motor neurons MN, and causing muscle fibre activation MF which accelerates shortening. Reflexes in distant segments i and i + j (|j| > 1) mutually inhibit one another via interneurons IN. This limits the number of moving segments to allow centre of mass motion (see text).

https://doi.org/10.1371/journal.pcbi.1006635.g002

Let us now consider the second criterion for peristaltic locomotion. Assuming all segment boundaries are of equal mass, the force on the centre of mass of the larva is proportional to the sum of the forces acting on each segment boundary, i.e. (2)

Newton’s third law tells us that any forces of interaction between segment boundaries (i.e. viscoelastic and muscle forces) must be of equal magnitude and opposite direction, so that they cancel in this summation and we are left only with contributions arising from substrate interaction. If the motion of the body is such that some number n_f of segments move forward at a given time, against a frictional force −μ_f, while n_b segments remain anchored or move backward, experiencing a frictional force μ_b, then the summation becomes (3)

In the limiting case of isotropic (direction-independent) substrate interaction we have μ_b = μ_f, and this expression tells us that the centre of mass will accelerate in the forward direction only when there are less segments moving forward than are moving backward or anchored to the substrate. Similarly, moving a small number of segments backward while the others remain anchored will result in backward acceleration of the centre of mass. Therefore, if the animal is to move relative to its substrate, it must ensure that only a limited number of its segments move in the overall direction of travel at a given time (indeed, this matches observations of the real larva [3, 22]). A more lengthy exposition of this requirement on limbless crawling behaviours can be found in [51].

We fulfil the requirement for a small number of moving segments by introducing mutually inhibitory interactions between the segmentally localised reflex circuits (Fig 2). We add a single inhibitory interneuron within each segment. When the sensory neuron within the local reflex activates, it excites this interneuron, which then strongly inhibits the motor neurons and inhibitory interneurons in non-adjacent segments, effectively turning off the local reflexes in distant neighbours. Adjacent segments do not inhibit each other in our model, allowing reflex activity to track mechanical disturbances as they propagate from one segment to the next. We comment on the plausibility of this feature of our model, given the experimental observation of nearest-neighbour inhibitory connections in the larval ventral nerve cord [28], in the discussion. Similarly, the head and tail segments do not inhibit each other, which permits peristaltic waves to be (mechanically) reinitiated at one extremity as they terminate at the other (see Discussion at the end of the previous subsection). This effectively introduces a ring-like topology into the neural model, matching our model of axial mechanics which couples head and tail motion through the total length constraint [40].

We now have a neuromuscular model consisting of four cell types repeated in each segment—sensory neurons, inhibitory interneurons, motor neurons, and muscle fibres. For the sake of simplicity we model all neurons as having a binary activation state governed by the algebraic relation (4) where V_i is the activation of the i’th cell, θ_i is its activation threshold, V_j is the activation of the j’th presynaptic cell, and w_j is the associated synaptic weight. Numerical values for the weights and thresholds used in our model are given in S1 Table, supplemental. Note that the muscle tension over a segment either vanishes (when the muscle fibre is in the inactive state) or has fixed magnitude b_i (when the muscle fibre is activated by local sensory feedback). For this reason we refer to b_i as the reflex gain.

Our choice to neglect neural dynamics is based on the large difference in timescales between the neural and mechanical dynamics. Typical motor neuron spiking occurs with a timescale on the order of 10⁻³ seconds. Spiking is observed to be significantly “averaged out” by the graded (non-spiking) muscle fibre responses, and respond on the order of ∼10⁻¹ seconds to prolonged motor neuron spiking [52, 53]. During locomotion, segmental compressions are driven by several longitudinal muscle fibres activating simultaneously [3, 14, 29] in response to largely independent motor neuron populations [54, 55] which are unlikely to spike with identical timing. This spatial integration should further “mask” the effects of neural dynamics. Note that the lack of neural dynamics in our model immediately rules out central pattern generation. However, this does not prevent our model from producing complex, larva-like behaviour, and we consider how our model could account for observations of central pattern generation in the discussion.

To summarise, the neural model we have constructed can be seen as consisting of two parts, a segmentally repeating local reflex and a mutual inhibition circuit acting between non-adjacent reflexes. The local reflex is constructed so that muscles will act as motors, amplifying segmental compressions and counteracting friction. The mutual inhibition circuit couples distant reflexes to allow only localised amplification. By limiting the number of moving segments, this should ensure that the model larva can produce a net force on its centre of mass.

Larva-like axial compression waves and lateral oscillations result from conservative mechanics

One of the advantages of grounding our model of larval exploration in the body’s physics is that we are now able to apply powerful analytical results from classical mechanics in order to understand the model’s behaviour. In this section we attempt to elucidate the naturally preferred motions of the larva by focusing our attention on the conservative mechanics of the body while neglecting friction forces, which would cause all motion to stop, and driving forces, which might impose a particular pattern of motion.

In this case, the general character of motion is specified by the Liouville-Arnold integrability theorem. This theorem asks us to look for a set of conserved quantities associated with a mechanical system, which remain unchanged as the system moves (energy, momentum, and angular momentum are examples of some commonly conserved quantities). If we can find a number of these quantities equal to the number of mechanical degrees of freedom in our model, then the theorem tells us that the motion of the body is integrable—it can be expressed analytically, and must be either periodic or quasiperiodic. If there are not enough conserved quantities, then the system is said to be nonintegrable. In this case the motion is much more complicated and will be chaotic for some initial conditions. These chaotic motions do not permit analytical expression and must be approximated through simulation.

In this section, we explicitly seek a case for which there is a “full set” of conserved quantities (we provide only major results here, for detailed derivations see S2 Appendix). We begin by restricting ourselves to considering only small deformations of the larval midline, in the case where all segments are of identical axial stiffness k_a, transverse stiffness k_t, mass m, and length l. Under these assumptions, the total mechanical energy of the body may be written (5) where x and y are vectors giving the displacement of each segment boundary along the body axis and perpendicular to the body axis, respectively, p_x and p_y give the translational momentum associated with each direction, D₂ and D₄ are difference matrices arising from a Taylor series expansion of our model’s potential energy (see S2 Appendix), and and are characteristic axial and transverse frequency scales. By making a linear change of coordinates {x, y, p_x, p_y} → {X, Y, p_X, p_Y} to the eigenbasis of D₂ and D₄ (see S2 Appendix) this simplifies to (6) where λ_a,i and λ_t,i are eigenvalues associated with the coordinate transformation. This expression is a sum of component energies, each of which is independently conserved. The Liouville-Arnold theorem immediately tells us that the motion of the body must be (quasi)periodic in the case of conservative small deformations. Indeed, the energy associated with each of the new coordinates X_i, Y_i is in the form of a harmonic oscillator, and each coordinate executes pure sinusoidal oscillations. By transforming back to the original coordinates x, y we obtain a set of collective motions (modes) of the body which are energetically isolated and have a sinusoidal time dependence, corresponding to axial and transverse standing waves. We will refer to the X_i, Y_i as modal coordinates since they describe the time dependence of each of the collective motions.

Each transverse standing wave corresponds to a periodic lateral oscillation of the body, with a unique frequency given by . We determined these frequencies numerically, along with the spatial components of the lowest frequency standing waves (Fig 3A). These can be seen to match the eigenmaggot shapes extracted from observations of unbiased larval behaviour [56].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Conservative, small-amplitude motions of the body decompose into a set of axial and transverse standing waves.

A: spatial component of first four transverse standing waves (top, black) compared to first four experimentally determined eigenmaggots [56] (top, blue), with natural frequencies of oscillation (bottom). B: spatial component of first four axial standing waves (top), with natural frequencies of oscillation (bottom). Note that axial standing waves come in pairs with identical frequency. C: Pairs of axial standing waves can be combined to produce forward-propagating (top, solid black line) and backward-propagating (bottom, solid black line) travelling waves. Head and tail extremities move in phase (dashed black line) due to our total length constraint (see text), reminiscent of the “visceral pistoning” observed in the real animal [3].

https://doi.org/10.1371/journal.pcbi.1006635.g003

The axial standing waves correspond to oscillating patterns of segmental compression and expansion. While each transverse standing wave had its own unique frequency of oscillation, the axial standing waves come in pairs with identical frequency but different spatial components—each member of the pair corresponds to a different spatial pattern of segmental compression/expansion spread across the body, but these patterns oscillate in time with the same frequency. We were able to analytically determine the frequency of oscillation of the i’th pair of axial standing waves to be (7)

The spatial components of the axial standing waves could also be determined analytically (8)

Where x_k,i is the displacement of the k’th segment boundary for the i’th pair of standing waves. We plot the axial frequencies ω_a,i and spatial components x_k,i in Fig 3B.

The fact that the axial oscillation frequencies come in identical pairs allows us to combine the axial standing waves with a ±90° relative phase shift to form pairs of forward and backward travelling wave solutions (see S2 Appendix for the full derivation) (9) where x_k,i(t) gives the displacement of the k’th segment boundary as a function of time for the i’th pair of travelling waves. The choice of a plus or minus sign corresponds to the choice between forward or backward wave propagation. These solutions correspond to propagating waves of segmental compression and expansion similar to those seen during larval peristalsis. We plot the lowest frequency pair of axial travelling wave solutions in Fig 3C, and directly visualise the synthesis of travelling wave solutions from standing wave solutions in S2 Video.

To summarise, in this section we have shown that for the case of conservative, small oscillations, the motion of the body may be decomposed into a combination of transverse standing waves and axial travelling waves. This is of clear relevance to understanding the behaviour of the larva, which moves across its substrate by means of axial peristaltic waves while reorienting using lateral oscillations. Indeed, the transverse modes of oscillation that we have derived here match principal components of bending computed from actual larval behaviour [56]. Our results can be interpreted as providing a physical basis for these observations—the principal components extracted from real larval data correspond to a “natural” coordinate basis that is grounded in the animal’s mechanics. Furthermore, the proportion of postural variance explained by each principal component of the experimental data decreases with increasing modal frequency in our model (and thus increasing energy). We can therefore help to explain the observed ordering of principal components, as this corresponds to the larva “preferring” to occupy low-frequency, low-energy modes during most of its behaviour. We comment further on this observation in S3 Appendix in the context of axial modes. We will now focus on the small-amplitude motion of the body in the presence of energy dissipation due to friction and driving forces.

Forward and backward peristaltic locomotion can be obtained from simple reflexes

Reintroducing friction will clearly lead the motions described above to eventually terminate due to energy dissipation, unless opposed by transfer of power. In a previous section (Models—Neuromuscular system, see also S1 Appendix), we introduced a neuromuscular system to produce power flow into the body, but as described, it can only directly transfer power into the axial degrees of freedom. In the small deformation model we have just analysed, the axial and transverse degrees of freedom are energetically decoupled. It follows that transverse friction is unopposed and any transverse motion must eventually terminate in the case of small deformations. In this section we will therefore focus only on the axial degrees of freedom, which correspond to the peristaltic locomotion of the larva.

In Fig 4, we show the effect of coupling our neuromuscular model to the axial mechanics. For small reflex gains, the power flow into the body from the musculature is too low to effectively counteract frictive losses and the larva tends towards its passive equilibrium state, in which there is complete absence of motion. However, increasing reflex gain past a certain critical value leads to the emergence of long-term behaviours in which the larva remains in motion, away from its passive equilibrium. These motions correspond to forward and backward locomotion, driven by forward and backward propagating compression waves (see below), as predicted from our earlier description of the conservative motions of the body, and as observed in the real larva [3]. The qualitative changes in behaviour that occur as reflex gain is varied are depicted in Fig 4A, where we have measured the long-term centre of mass momentum of the larva, along with the long-term relative phase of the lowest frequency modal coordinates.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Emergence of limit cycles for forward and backward locomotion in the dissipative, small-amplitude model.

A: as reflex gain is increased past a critical point, the model larva attains a positive or negative long-term average centre of mass momentum (top, red and blue lines), signifying continuous forward or backward motion relative to the substrate, and a ±π/2 relative phase difference between the two lowest frequency axial standing wave modes (bottom, red and blue lines), signifying the presence of forward- or backward-propagating axial travelling waves. B: trajectories of individual point masses in the model for forward (top) or backward (bottom) locomotion (see S2 Fig for corresponding neural state). C: projection of model trajectories onto the lowest frequency axial modes and the centre of mass momentum reveals a pair of (putative) stable limit cycles for forward (blue) and backward (red) locomotion. Parameters used to generate this figure are given in S2 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g004

The exact value of reflex gain required for onset of locomotion depends on the particular mechanical parameters used in our model (see S2 Table for parameters used in Fig 4). In principle, this bifurcation point of the dynamics should be amenable to analytical investigation by studying the stability of the linearised model dynamics around the passive equilibrium state [57, 58]. In practise, however, the presence of hard nonlinearities in the sliding friction model makes such an approximation inaccurate.

For low reflex gains the centre of mass momentum tends to 0 as the body comes to a stop and enters a passive equilibrium state. The relative phase of the low frequency modal coordinates tends to either 0 or 180 degrees, which also corresponds to a loss of momentum. For larger values of reflex gain, the total momentum is either positive, zero, or negative. Positive and negative total momentum represent forward and backward locomotion, respectively, while zero momentum corresponds to two unstable cases which we discuss below. The relative phase of the lowest frequency modal coordinates tends to ±90° at high reflex gains, corresponding to the presence of forward- or backward-propagating compression waves (see previous section). As in the larva [1, 3], forward-propagating waves drive forward locomotion while backward-propagating waves drive backward locomotion (Fig 4B).

We believe that these behaviours arise from the production of a pair of limit cycle attractors in the system’s phase space, which we visualise in Fig 4C by projecting the system state onto the lowest frequency pair of axial modes, and plotting the associated modal coordinates along with the centre of mass momentum. Since wave motion implies that pairs of modal coordinates should perform pure sinusoidal oscillations with equal amplitude and frequency, and a ±90° relative phase shift (see previous section and S2 Appendix), the travelling wave trajectories of the system become circles in this coordinate system (see discussion of Lissajous figures, [59]). Forward and backward locomotion can then be distinguished by the momentum of the centre of mass.

In this model, the speed of forward and backward locomotion are equal for a fixed value of reflex gain, while in the real larva the speeds are known to differ [15]. We comment on some possible explanations for this difference in the discussion.

In S2 Fig we show the neural state of the model larva during forwards and backwards locomotion. As expected given our previous exposition, we observe waves of activity in the nervous system which track the mechanical waves propagating through the body. Our sensory neurons also show a second, brief period of activation following propagation of the mechanical wave caused by a slight compression that occurs as segments return to equilibrium. This activity is “cancelled out” by the mutual inhibition circuit, so that motor neurons do not exhibit a secondary burst of activity.

The larva has zero long-term total momentum in the presence of large reflex gain in only two cases, both of which are highly unstable. First, if we start the larva so that it is already in its passive equilibrium state, so that no relative motion of segment boundaries occurs, it is obvious that there will be no activation of local reflexes and the larva will not spontaneously move out of equilibrium. The second case corresponds to a pure axial standing wave. In this case the larva is divided into two regions by nodal points where the axial displacement is zero. During the first half-cycle of the standing wave, one region will experience compression while the other experiences expansion, and in the second half-cycle these roles will reverse. The neural circuit we have constructed can amplify compression during both half-cycles since they are separated by a configuration in which no compression or expansion occurs, and this allows the entire nervous system to become inactive and “reset”. Since these behaviours are extremely unstable and require very specific initial conditions to be observed, we have not visualised them here.

While the mutually inhibitory connections in our model are not required for the propagation of locomotor waves, which can be maintained entirely by local reflexes [40], these connections do greatly enhance stability. In the absence of the mutual inhibition circuit, small mechanical disturbances in any stationary body segments can be amplified, giving rise to multiple compression waves which travel through the body simultaneously. If this instability produces an equal number of forward and backward moving segments then overall motion of the larva relative to the substrate will stop, in line with the argument presented earlier. We have also observed that roughly symmetrical substrate interaction is required to produce both forward and backward locomotion in our model. If friction is too strongly anisotropic, then locomotion can only occur in one direction regardless of the direction of wave propagation.

It is worth noting that the axial model presented in this section does display discrete behavioural states. However, there are no explicit decisions regarding which behavioural states to enter, since the particular neural states occupied during forwards and backwards locomotion are indistinguishable. Rather, both the apparent decision and the eventual direction of travel can only be understood by examining the dynamics and mechanical state of the body.

Conservative chaos from mechanical coupling

Having successfully produced peristaltic locomotion using our model, we will now turn our attention to the larva’s turning behaviours. As before, we will start from physical principles. In a previous section (Results—Conservative axial compression waves and transverse oscillations) we showed that, for the case of conservative small oscillations, transverse motions of the body were energetically decoupled from axial motions, and could be decomposed into a set of periodic standing waves. We will first extend our previous analysis to the case of energy-conservative, large amplitude motions in the absence of damping and driving; and then in the following section consider the impact of energy dissipation and the neuromuscular system on transverse motion.

To keep our presentation simple and allow visualisation of model trajectories, we will focus on a reduced number of the mechanical degrees of freedom. Namely, we will examine the bending angle ϕ and axial stretch q of the head segment (Fig 5A). We introduce an amplitude parameter ϵ by making the substitutions q → ϵq and ϕ → ϵϕ, so that the total mechanical energy of the head may be written in nondimensional form as (see S4 Appendix) (10) where p_q, p_ϕ are the radial and angular momentum of the head mass, and we have scaled the time axis of the model so that the natural frequency of axial oscillation is unity. The parameter λ = ω_t/ω_a = k_t/k_al² is the ratio of transverse and axial frequencies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. A reduced model of large amplitude motion.

A: we focus on the conservative dynamics of the head’s strain q and bend ϕ coordinates as amplitude ϵ is varied. B: head trajectories are visualised by Poincare section, in which the head’s configuration q, ϕ is plotted at successive turning points of the transverse bending motion (at which angular velocity vanishes, dϕ/dt = 0).

https://doi.org/10.1371/journal.pcbi.1006635.g005

In the case of small oscillations, i.e. ϵ → 0, the mechanical energy reduces to the simpler expression (11) which is clearly a sum of independent axial and transverse energies. These energies are individually conserved, so that the Liouville-Arnold theorem applies, and the motion of the head is integrable and (quasi)periodic. This is in clear agreement with our earlier results. For the more general case of large amplitude motion (ϵ > 0), we may have in principle only a single conserved quantity—the total energy of the system. Indeed, it should be clear from the presence of a “mixed” term in the mechanical energy (Eq 10) that the axial and transverse motions are energetically coupled at large amplitudes, so that the individual energies associated with each motion are no longer independently conserved. Given that we have only one conserved quantity for a two degree of freedom system, we can no longer rely on the Liouville-Arnold theorem to prove (quasi)periodicity of the motion, and must accept that the system’s behaviour may be chaotic.

To investigate this possibility further, we first note that conservation of energy implies that flow within the four dimensional phase space must be constrained to lie on the energy surface given implicitly by the relation H(q, ϕ, p_q, p_ϕ) = E. Therefore, given a particular value E for the total energy, the system dynamics becomes three dimensional. This allows us to visualise the behaviour of the system by plotting the points at which trajectories intersect a two-dimensional Poincare section [57, 58]. We define our Poincare section by the condition that the angular momentum vanishes p_ϕ = 0 (equivalently, angular velocity vanishes dϕ/dt = 0), and plot successive crossings of the section as points in the q, ϕ plane. This allows us to intuitively interpret points in the Poincare section as configurations of the head at successive turning points (extrema) in the transverse motion (Fig 5B).

In what follows, we set the total energy to be so that when ϵ = 1 we can in principle obtain complete compression of the head (q = −1). We choose to set , giving an irrational frequency ratio. This loosely matches observations of the real larva in which the frequency of transverse oscillations is approximately half that of axial oscillations but does not satisfy an exact (rational) resonance relationship [44]. The results we obtain with these parameters do not differ much from results for other energies or other frequency ratios, including resonant relationships.

Poincare plots for the cases ϵ → 0 and are shown in Fig 6. When ϵ → 0 (Fig 6A), conservation of transverse energy implies that the turning points of the transverse motion must remain constant. The fact that the frequency ratio λ is irrational implies that the overall motion is quasiperiodic, and the values of q obtained at successive transverse turning points should not repeat. In accordance with these observations, the Poincare section for ϵ → 0 consists of a series of verticle lines (Fig 6Ai). For the Poincare plot becomes distorted, but the majority of trajectories still trace out one-dimensional curves in the section (Fig 6Bi), which is indicative of persistent quasiperiodic behaviour. At the Poincare plot appears qualitatively different. There is now a large region of what appears to be “noise”, indicating that the configuration of the head at successive transverse turning points has become unpredictable. This is a clear signature of deterministic chaos. There are, however, several regions of the section indicative of (quasi)periodic behaviour. These appear as one-dimensional curves or discrete points in the Poincare section (Fig 7Ai). At ϵ = 1, the region of the Poincare plot occupied by chaos has expanded, although there still appear to be some regions corresponding to (quasi)periodic behaviour (Fig 7Bi). These results qualitatively agree with the Kolmogorov-Arnold-Moser theorem [59], which tells us that quasiperiodic behaviour should persist under small nonintegrable (chaotic) perturbations of an integrable Hamiltonian, and that the region of phase space corresponding to chaotic behaviour should grow with the perturbation size (in our case, the perturbation size corresponds to the amplitude of motion ϵ). However, our model as presented here does not formally meet the requirements of this theorem (see S4 Appendix).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Emergence of deterministic chaos in the conservative head dynamics as amplitude of motion is increased.

A, B: for small amplitudes (ϵ → 0, ϵ = 1/3), Poincare section shows quasiperiodic head oscillations (i), while the maximum Lyapunov characteristic exponent (MLCE), which quantifies the dominant rate of separation of nearby phase trajectories, converges to ∼ 0bits s⁻¹ (ii), the power spectra of head stretch q and bend ϕ coordinates show clear peaks with little “noise” component (iii), and autocorrelations of these variables decay linearly (iv). These results betray non-chaotic, quasiperiodic oscillations for small amplitudes. MLCE, power spectra, and autocorrelations were computed for initial conditions shown by black dot in panel i. Parameters used to generate this figure are detailed in the main text, and reported in S3 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. (continued from Fig 6) Emergence of deterministic chaos in the conservative head dynamics as amplitude of motion is increased.

A, B: for large amplitudes (ϵ = 2/3, ϵ = 1), the Poincare section contains a large chaotic sea (i), while the MLCE converges to a positive value (ii), power spectra become “noisy” (iii), and autocorrelations decay rapidly (iv). These results strongly suggest the existence of deterministic chaotic head dynamics for large amplitudes. MLCE, power spectra, and autocorrelations were computed for initial conditions shown by black dot in panel i. Parameters used to generate this figure are detailed in the main text, and reported in S3 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g007

Analysis by Poincare section provides an invaluable method to determine the character of overall system behaviour by direct visualisation of trajectories, but cannot be applied to higher dimensional systems. This is problematic, since we would like to assess the existence of chaos beyond our reduced model of the larva’s head. We therefore deployed a series of other methods which are possibly less reliable than the method of Poincare section but can be applied equally well to higher dimensional systems. These included estimation of the maximal Lyapunov characteristic exponent (MLCE) for the dynamics along with calculation of the power spectrum and autocorrelation of internal variables [57, 58, 60]. The MLCE can be thought of as quantifying the rate of separation of nearby trajectories, or, equivalently, the rate at which information is generated by the system being analysed [61]. A positive MLCE is generally considered a good indicator of chaotic behaviour. The power spectrum of a periodic or quasiperiodic process should consist of a “clean” set of discriminable peaks, whereas that of a chaotic process should appear “noisy” and contain power across a wide range of frequencies. Meanwhile, the autocorrelation of a periodic or quasiperiodic process should show a strong oscillatory component with an envelope that decays linearly with time, while that of a chaotic process should show a much quicker decay, similar to a coloured noise process. In Fig 6 we plot these measures at each value of ϵ, for a trajectory starting with initial conditions indicated on the corresponding Poincare plot by a filled black circle. These measures confirm increasingly chaotic behaviour as the amplitude ϵ increases, in agreement with our Poincare analysis. In Fig 8 we show a solution including all degrees of freedom in our conservative mechanical model (i.e. not just those of the head). The trajectory of individual segments relative to the substrate appears qualitatively irregular, while the indirect measures we introduced above (MLCE, power spectrum, autocorrelation) are all indicative of deterministic chaotic behaviour.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Conservative planar motion of the body is chaotic at large amplitudes.

A: trajectories of individual segment boundaries appear qualitatively irregular, B: our estimate of the maximum Lyapunov characteristic exponent converges to a positive value, C: power spectra of head stretch q and bend ϕ show a strong “noise” component, and D: their autocorrelations decay rapidly. All are indicators of deterministic chaos. Parameters used to generate this figure are given in S4 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g008

Spontaneous turning and reversals require no additional control

We will now reintroduce dissipative effects into our model of larval motion in the plane. We previously saw that conservative mechanics predicted axial travelling waves and transverse oscillations. These were lost in the presence of friction, but the axial travelling waves could be recovered with the addition of a neuromuscular system designed to selectively counteract frictive effects. We have now seen that conservative mechanics predicts chaotic planar motion. Although our neuromuscular model transfers power only into the axial degrees of freedom, we recall from the previous section that axial and transverse motions are energetically coupled at large amplitudes. We therefore tentatively expect that we may be able to recover the complete chaotic planar motion without requiring any additional mechanism for direct neuromuscular power transfer into tranverse motion.

To investigate whether our dissipative planar model shows chaotic behaviour, we ran n = 1000 simulations starting from almost identical initial conditions (euclidean distance between initial mechanical state vectors < 10⁻⁷, with no initial neural activity) and identical parameters (see S5 Table). We can indeed observe that the simulated larva perform peristalsis with interspersed bending of the body (turns), and that the fully deterministic system produces apparently random turning such that the simulations rapidly diverge (S3 Video). Since most working definitions of chaos require strictly bounded dynamics, we here restrict our analysis to the coordinates describing deformation of the body (segmental stretches and bending angles), ignoring motions of, or overall rotations about, the centre of mass (i.e., the trajectory through space of the body, which we will analyse in the following section).

Qualitatively, the deformations of the large amplitude dissipative model appear irregular (Fig 9A). However, there are persistent features reminiscent of the ordered small-amplitude behaviour described in previous sections. In particular, there are clear axial travelling waves and lateral oscillations. However, whereas forward- and backward-propagating axial waves previously corresponded to stable limit behaviours, the large amplitude system appears to go through occasional “transitions” between these behaviours. In addition, apparently spontaneous large bends appear occasionally, suggesting that the neuromuscular system is effectively driving transverse motion.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Dissipative planar motion is chaotic.

A: representative segmental stretch (left) and bend (right) time series (see S3 Fig for corresponding neural state). Note the occurence of a large bend starting at ∼ 1–2 seconds at the larva’s head, which appears to propagate backwards along the body while triggering a “transition” from forward to backward wave propagation at ∼ 3.5 seconds. Forward wave propagation resumes at ∼ 6 seconds. B: power spectra of the head stretch q (top) and bend ϕ (bottom) showing significant “noise” component. C: Autocorrelations of q and ϕ rapidly decay. D: probability density of correlation dimension estimates for 1000 mechanical trajectories. The dimension of the system’s limit set is estimated as ∼ 3.5 (median, vertical blue line). E: maximum Lyapunov characteristic exponent estimates converge to a positive value. All measures suggest the presence of deterministic chaotic dynamics. Parameters used to generate this figure are given in S5 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g009

The irregularity of the axial motion is reflected in the pattern of sensory neuron activation (S3 Fig). However, the mutual inhibitory interactions in our model act to filter this input, allowing only a small window of excitability within the central nervous system. As a result, interneuron and motor neuron activity appears fairly ordered, tracking and reinforcing axial compression waves.

We used four measures to assess whether our qualitative observation of irregular motion betrays the existence of deterministic chaos. First, we analysed the power spectrum of individual cooordinates (Fig 9B). The power spectra of all degrees of freedom showed consistent harmonic peaks along with a strong “noisy” non-harmonic component, a pattern consistent with chaos and incommensurate with (quasi)periodicity (Fig 9B shows data for head bending angle and stretch Next, we computed the autocorrelation of the same degrees of freedom. The autocorrelations of all degrees of freedom showed a periodic pattern with a peak at 0 seconds time lag followed by a rapid decay (Fig 9C). This is characteristic of oscillatory chaotic behaviour, and the exponential loss of information regarding initial conditions that chaotic systems display. We then estimated the correlation dimension (D_c) of the limit set of our model’s dynamics. Note that we did not apply this measure to the conservative models in the previous section since the Liouville theorem rules out attracting limit sets for these systems. The distribution of correlation dimension estimates for our dissipative system across all 1000 trials is shown in Fig 9D. Estimates were clustered around ∼ 3.5 (median dimension 3.46), with 93% of estimates lying in the range [3–4]. These results are suggestive of a limit set with fractal dimension, which is a signature of “strange” chaotic attractors. Finally, we computed an estimate of the maximal Lyapunov characteristic exponent (MLCE). As in the previous section, we estimated the MLCE for our system to be positive (∼ 13textrmbitss⁻¹, Fig 9D), a very strong indicator of chaotic behaviour. All of these results point to the presence of a chaotic dynamical regime in our dissipative large amplitude model.

Exploration emerges without decisions or stochasticity

As the coupled biomechanical and neuromuscular system produces both forward and backward peristalsis and lateral oscillations, each simulated larva will trace out a 2D trajectory over time. As shown in Fig 10A, the chaotic deformations characterised in the previous section caused the larvae to disperse across their two-dimensional substrate, following paths reminiscent of the real animal’s exploratory behaviour.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Deterministic exploration.

A: dispersion of the centres of mass of 1000 simulated larvae, starting from almost identical mechanical initial conditions (overlayed in inset). B: tortuosity and fractal (box-counting) dimension for all 1000 paths indicate plane-filling behaviour (blue line = mean tortuosity, red line = mean dimension, see text; see S4 Fig for power law analysis of trajectory curvature and angular speed). C: mean-squared displacement (black line) shows transient quadratic growth (blue line) followed by asymptotic linear growth (red line, asymptotic diffusion constant ≈ 144segs²s⁻¹; see also log-log plot, S5 Fig). D: distribution of body bends (black) with maximum likelihood von Mises (red) and wrapped Cauchy (blue) fits. E: distribution of run lengths with maximum likelihood exponential fit (red). Run lengths were calculated as duration between successive crossings of a threshold body bend (20°), indicated by blue lines in panel D. See S6 Fig for analysis of tail speed and head angular velocity. Parameters used to generate this figure are given in S5 Table.

https://doi.org/10.1371/journal.pcbi.1006635.g010

To characterise the trajectories of our model, we first investigated them at a global level, based on the centre of mass (COM) trajectory of each simulated larva, computing the tortuosity and fractal dimension of the paths (Fig 10B) [62]. We defined our tortuosity measure as (12) where D is the net displacement of the COM between initial and final times, and L is the total distance travelled by the COM along its path. Note that if the COM travels in a straight line between initial and final times we will have D = L so that T = 0. In the limit L → ∞, corresponding to the COM taking an extremely long path between its initial and final states, we have so that T → 1. We calculated the mean tortuosity of our COM trajectories to be 〈T〉 = 0.43, with a variance of 〈(T − 〈T〉〉)² = 0.05. The lowest (highest) tortuosity observed was T = 0.05 (T = 0.95).

We estimated the fractal dimension D_b of the COM trajectories using a box-counting algorithm. The minimum expected dimension D_b = 1 would correspond to curvilinear paths (e.g. straight line or circular paths), while the maximum expected dimension of D_b = 2 corresponds to plane-filling paths (e.g. brownian motion). We calculated the mean dimension of the COM trajectories to be 〈D_b〉 = 1.37, with variance 〈(D_b − 〈D_b〉)²〉 = 0.01. The lowest (highest) path dimension observed was D_b = 1.17 (D_b = 1.95). We have plotted the tortuosity and fractal dimension of every path in Fig 10B. These results show that the trajectories of the model differed markedly from straight lines (tortuosity T > 0 and box-counting dimension D_b > 1), and displayed a good ability to cover the planar substrate (box-counting dimension 1 < D_b < 2). We also note that our COM trajectories display the power-law relationship between angular speed and curvature reported by [63], with a scaling exponent (β ≈ 0.8) falling within the range reported for freely exploring larvae (S4 Fig).

We next investigated the rate at which the simulated larvae explored their environment. To do this, we calculated the mean-squared displacement (MSD) of the COM over time (Fig 10C). This is a standard measure used to characterise diffusion processes, and is defined as (13) where R_i(t) is the position of the i’th larva’s COM at time t and n = 1000 is the number of trials being averaged over. We observed an initial transient, lasting on the order of 10 seconds, during which the MSD grew as ∼ t² (blue line, S5 Fig), after which growth slowed and tended to ∼ t (linear fit for t > 80 seconds shown by red line, Fig 10C and S5 Fig, r² = 0.99, diffusion constant D = 144segs²s⁻¹). The initial transient was not due to our particular initial conditions, since it remained even after discarding > 60s of initial data. These results therefore tell us that, although on long timescales our model appears to execute standard Fick diffusion or a Brownian random walk (linear growth of MSD), on short timescales the model’s behaviour is superdiffusive (approximately quadratic growth of MSD) [64, 65]. This is in good agreement with observations of the real larva [6, 48]. The superdiffusive behaviour of the larva was previously explained in terms of a persistent random walk [6], in which the larva’s current and previous headings are highly correlated during straight runs so that the animal follows an approximately ballistic trajectory on short timescales. We believe that persistence effects arise in our model due to the finite time required for the deterministic chaotic dynamics to destroy information regarding initial conditions.

We next calculated some other standard measures found in the larva literature. To do so, we built a two-segment representation of each simulated larva by drawing vectors from the tail extremity to the A2-A3 segment boundary (the tail vector, T), and from the A2-A3 boundary to the head extremity (the head vector, H). We then defined the body bend, θ, to be the angle between tail and head vectors, θ = atan(H_y/H_x) − atan(T_y/T_x). The head angular velocity ν was computed as the cross-product of the head vector and the head extremity’s translational velocity measured relative to the lab frame, , while the tail speed v was taken to be the magnitude of the tail extremity’s translational velocity measured relative to the lab frame, . The tail speed and head angular velocity both show a strong oscillatory component, which can be seen in the time and frequency domains (S6 Fig). The power spectra of these variables contains considerable “noise” over a wide spread of frequencies, in accordance with the results of the previous section. The distribution of tail speeds for our model is bimodal, similar to that of the real larva [44]. The body bend angle was observed to be symmetrically distributed (Fig 10D), with roughly zero mean (〈θ〉 = 0.005), small variance (〈(θ − 〈θ〉)²〉 = 0.13), slight positive skew (SK(θ) = 0.23), and high excess kurtosis KU(θ) = 7.9. The kurtosis of our data precludes a good fit by the von Mises distribution (maximum likelihood estimate shown by red line in Fig 10D). Our data appears to be better fitted by a wrapped Cauchy distribution, though this overestimates the central tendency of our data (maximum likelihood estimate shown by blue line in Fig 10D). The high excess kurtosis of the body bend distribution gives a quantitative measure of the rare large bending events mentioned at the beginning of the previous section, and qualitatively matches experimentally observed distributions of real larval bends [44, 66]. Our model also reproduces the observed overall speed of larval locomotion (median model speed = 0.26 body lengths s⁻¹, real larval range ∼ 0.1 − −0.35 body lengths s⁻¹), the turn rate (median model turn rate = 2.08 turns min⁻¹, real larval range ∼ 0 − −4.5 turn min⁻¹, threshold body bend for turn classification = 30^deg to match relevant literature), and the relative distance gained during free locomotion (median model distance gained = 0.14 body lengths s⁻¹, real larval range ∼ 0 − −0.2), with our results being more consistent with observations of third instar than first instar larvae [66].

Finally, we computed a run-length distribution by setting a threshold body bend angle θ_turn = 20° (as in [13]) and calculating the length of time between successive crossings of this threshold. The distribution we obtained appears approximately linear on a log-linear plot (Fig 10E, linear fit r² = 0.99 with slope λ = −0.075), and is better fit by an exponential than a power law distribution (maximum likelihood estimates, log likelihood ratio = 5281, p < 0.01). Together with our observation of asymptotic linear growth of MSD, the exponential distribution of run lengths suggests that the model can be considered to be effectively memoryless on long timescales [65]. This again agrees with the observed rapid loss of information from the system due to its chaotic dynamics, as quantified by the Lyapunov exponent and autocorrelation analysis of the previous section.

Ultimately, the analysis of our model supports a view of the larval exploratory routine as an (anomalous) diffusion process arising from the deterministic chaotic dynamics of the body. The model nervous system functions purely to recover these dynamics from the effects of frictive energy dissipation, and to ensure centre of mass motion, rather than explicitly directing exploration.