Introduction

Drosophila larvae express a fairly tractable behavioral repertoire that is consistent across the 4-5 days of the larval life stage (Almeida-Carvalho et al., 2017). Most of the larval lifetime is dedicated to foraging for suitable nutrients while avoiding threats via stereotyped evasive behaviors. Behavioral control is achieved via neural circuits conserved throughout development (Gerhard et al., 2017), that span the entire tripartite CNS consisting of the brain, the subesophageal zone (SEZ) and the ventral nerve cord (VNC), making the larva a formidable system for studying control, execution and adaptation of behavior (Zarin et al., 2019; Jovanic, 2020; Eschbach and Zlatic, 2020; Imambocus et al., 2022). After reaching critical mass for pupation, homeostatic signals switch behavior towards food aversion, hypermobility and collaborative burrowing (Wu et al., 2003), terminating the feeding state and leading to pupation and metamorphosis.

High resolution methods for behavioral tracking (Ohyama et al., 2013; Risse et al., 2017; Schumann and Triphan, 2020; Tadres and Louis, 2020; Thane et al., 2023; Laurent et al., 2024), now routinely used in neuroethological experiments, have revealed detailed insight in the organisation of larval foraging behavior. In the absence of food resources, larvae engage in free exploration to locate food patches (Sims et al., 2019) with a characteristic alternation of locomotory activity and brief pauses (Sakagiannis et al., 2020), a property also reported for adult fly behavior (Ueno et al., 2012; Reynolds et al., 2015). Food consumption involves repetitive feeding motion and digging into the substrate (Kim et al., 2017). Statistical regularities that govern foraging behavior have been un-veiled by analysis at the microscale of body kinematics and at the macroscale of larva trajectories (Denisov et al., 2013; Risse et al., 2017; Karagyozov et al., 2018). Locomotion combines the basic sensorimotor primitives of crawling and turning, and has been in the main focus of recent studies (Heckscher et al., 2012; Wystrach et al., 2016; Sims et al., 2019), whereas studies of feeding behavior remain scarce (Ruiz-Dubreuil et al., 1996). Both, crawling and feeding consist of recurring sensorimotor cycles controlled by central pattern generating circuits (Heckscher et al., 2012; Mantziaris et al., 2020; Miroschnikow et al., 2018a).

Salient olfactory cues can trigger appetitive or aversive chemotaxis, during which larvae navigate up or down a chemical gradient (Gomez-Marin et al., 2011; Slater et al., 2015; Schleyer et al., 2015a; Klein et al., 2017). During appetitive chemotaxis, the detection of minor concentration changes during lateral bending causes a directional bias in the turning movement, establishing a mechanism of active sensing (Gomez-Marin and Louis, 2012; Wystrach et al., 2016; Thane et al., 2019). Encounters with novel odorants in the presence of a food reward or a negative reinforcement such as the bitter taste substance quinine can induce olfactory learning, enabling short- and long-term behavioral adaptations (Schleyer et al., 2011; Schleyer et al., 2015b; Gerber and Stocker, 2007; Diegelmann et al., 2013; Widmann et al., 2018; Weiglein et al., 2019; Jürgensen et al., 2024). For quantification of choice behavior, widely-used standard group assays for behavioral preference testing have been established (Gerber and Stocker, 2007; Gerber et al., 2014; Schleyer et al., 2015a; Schleyer et al., 2015b; Schleyer et al., 2011).

The routine access to detailed behavioral data, and the broad approaches to neuroethological experiments makes the Drosophila larva an ideal model system for computational modeling studies on behavioral control principles. Existing generative models typically aim at specific aspects of larval behavior and can vary widely in model type and abstraction level, ranging from basic neuromechanics to the abstract mathematical description. We therefore propose the concept of the behavioral architecture, a three-layered hierarchical and modular framework, nested from behavioral primitives to high-level function, to integrate models of diverse types and different levels of abstraction (Figure 1A). At the basic layer we establish a refined coupled-oscillator model for larval locomotion capable of autonomous exploratory behavior. To this end we provide a detailed analysis of larval tracking data from experiments of free exploration. Introducing sensory modulation at the second layer allows us to reproduce several experimental findings of appetitive chemotaxis. The combination with a biologically realistic spiking neural network model for associative learning at the adaptive layer reproduces empirical population-level learning scores across a multiple trial training protocol. The architecture allows for agent-based, potentially closed-loop, simulations of virtual larvae in virtual environments and generates simulated data at the level of the experimentally observed behavior facilitating model calibration and evaluation.

Figure 1.

Model

Behavioral architecture

Larval behavior is hierarchically structured, i.e. sensorimotor primitives such as crawling, bending and feeding can be integrated into more complex behaviors such as exploration, taxis or exploitation. It has been proposed that the hierarchy of animal behavior is reflected by the underlying neuroanatomy as a hierarchy of nested sensorimotor loops (Prescott et al., 1999; Wilson and Prescott, 2022; Prescott and Wilson, 2023). A functional modeling paradigm that exploits this idea regards the neural system as a layered control architecture (or subsumption architecture) (Brooks, 1986) where low-level stereotyped reflexive and repetitive behaviors are autonomously generated by localized peripheral motor circuitries, while multisynaptic loops involving more centralized neural circuits act as descending modulators on the local circuits in order to coordinate global and complex behavior. The central idea is that energy-efficient decentralized neural control is the rule, while higher centers are recruited only when more extensive integration is needed e.g. in order to react suitably to unexpected sensory stimulation. Furthermore, there are only limited degrees of freedom by which higher layers can influence local sensorimotor loops via descending pathways e.g. by starting or stopping, or by accelerating, decelerating or inverting their autonomous function (Sen et al., 2017; Feng et al., 2020; Bidaye et al., 2020). The complexity of inter-layer top-down modulation is thus predicted to be considerably lower than the complexity of signal integration within layers. Layered control architectures have been used mainly in behavior-based robotics (Bicho, 1999), allowing sequential calibration and modular integration of partial neuroscientific models under a common framework (Brooks, 1986; Prescott et al., 1999).

We here propose a three-layered behavioral architecture for Drosophila larva foraging as illustrated in Figure 1A. The bottom layer consists of three basic behaviors: crawling, turning and feeding. Each is realized by a low-level local sensorimotor loop between motor effectors and sensory feedback. For exploration this involves proprioception and mechanoreception, for feeding it additionally involves gustatory input. Integration of these basic behaviors within the layer gives rise to composite behaviors. For example while exploration in stimulus-free conditions only entails crawling and turning, differences in the temporal microstructure between these two yield a spectrum from localized search to remote dispersal. The peripheral ventral nerve cord (VNC) circuits mediating exploration have been shown to be capable of autonomous, decentralised function without the need for any top-down regulation (Sims et al., 2019).

The intermediate layer introduces salient sensory stimulation of different modalities and therefore allows for reactive behavior in the face of presented risks and opportunities. Following the subsumption architecture paradigm, we only assume a single modulatory input from the intermediate to the bottom layer. This assumption is in line with the idea of pre-motor integration of signals from different sensory modalities into a final integrated descending modulatory pathway, which eventually affects behavior (Wystrach et al., 2016; Eschbach and Zlatic, 2020). Modulation of exploratory behavior in the presence of salient stimulation enables coherent navigation along sensory gradients. In the present study we consider odor-driven chemotaxis as a process of active sensing in which the larva constantly samples the odor concentration during lateral bending motions, enabling odorscape navigation and odor source discovery. In an alternative application, we recently simulated larval behavior during active sensing in thermotaxis experiments across different Drosophila species (Kafle et al., 2024).

Finally, at the top layer, associative learning biases the evaluation of recurring sensory stimuli and can modulate the innate valence of odors to allow for experience-dependent adaptation of navigation. Again, the modulatory top-down signal to the underlying layers is limited to the modulation of approach versus avoidance behavior, in line with the research on the mushroom body output signal in Drosophila (Gerber and Stocker, 2007; Gerber et al., 2014; Schleyer et al., 2015a; Schleyer et al., 2015b; Schleyer et al., 2011; Owald and Waddell, 2015).

The proposed behavioral architecture is naive to the level of model abstraction of the underlying neural mechanisms. Its role is to create a modular framework where diverse and potentially competing models can be positioned, integrated and tested in behavioral simulations. It is therefore justified to use simple generative models to bridge the gaps that have not yet been studied in detail until they can be substituted by more elaborate future versions. To ensure flexibility and compatibility, every module has a defined basic input and output regardless of the internal structure, complexity, level of abstraction and spatiotemporal scale of its possible instantiations. For the purpose of modeling under the behavioral architecture framework we provide the larvaworld python package available at https://pypi.org/project/larvaworld/.

Locomotory model

We propose a model for locomotion in the two-dimensional plane that defines the momentary body state by the instantaneous forward velocity v and angular velocity ω as generated through crawling and bending, respectively. Extending on previous models (see Discussion) we incorporate two new features based on our kinematic analysis of larval locomotion: (i) the peristaltic cycle phase dependent attenuation of angular motion, and (ii) the intermittent crawling as transitions between runs and pauses. We briefly describe the modular components of our model. A detailed mathematical description is provided in the Materials and Methods section. The complete pipeline for model calibration is described in Appendix 2.

We choose to simplify the larva to a two-segment body (Figure 2F). This abstraction allows to describe the body state at any point in time by only three variables: (i) position of a selected spine midpoint (Appendix 1-Figure 1), (ii) absolute orientation of the front segment θ, and (iii) bending angle θ_b between the front and rear segments. This approach is in line with the common practice of quantifying larva bending via a single angle between one anterior and one posterior vector (Gomez-Marin et al., 2011; Lahiri et al., 2011; Paisios et al., 2017).

Figure 2.

Locomotion is achieved via forward crawling and lateral body bending. Crawling moves the midpoint along the orientation vector. Bending reorients the front segment by rotation around the midpoint. Forward displacement restores θ_b back to zero by gradually aligning the rear segment to the orientation axis. A demonstration of the bisegmental model simplification is shown in Video 4.

Based on our empirical results (Figure 2 and Figure 3), we propose a coupled oscillator model of locomotion that generates both, weathervaning and headcasts. To this end we implement an oscillatory process (crawler) that generates the forward-velocity v during subsequent stride cycles (Figure 1B) where a tonic input I_C modulates crawling frequency f_C . A second oscillator (turner) generates alternating left-right bending as previously proposed by (Wystrach et al., 2016) with oscillatory amplitude A_T . To incorporate the empirically observed crawl-bend interference, we couple the two oscillators by imposing a crawler-phase dependent suppression on A_T where the phase-dependent modulation is modelled as Gaussian function and fitted to the empirical data (red profile in Figure 3C).

Figure 3.

To accommodate the intermittent behavior with alternating crawling runs and crawling pauses, we chose to implement a stochastic model approach where the number of strides per stridechain and the pause durations each follows a log-normal distribution fitted to the empirical data (Figure 3B). For a demonstration of the locomotory model in different configurations see Video 1.

Under the constraints of the subsumption architecture paradigm the crawler and turner can be influenced by higher-order circuits in a limited number of ways. In the case of the crawler, frequency modulation and the initiation or cessation of the oscillation-cycle are the two available top-down modulations. Likewise for the turner, top-down modulation of the tonic input affects both, the frequency and the amplitude of oscillation.

Results

Kinematic analysis of larva locomotion

We start out with the kinematic analysis of experimental larva trajectories and body postures in order to infer and parameterize several aspects of larva locomotion that will inform our modeling approach. Using diverse metrics that capture spatial and temporal dynamics we specifically assess the oscillation of forward velocity during consecutive peristaltic strides and its influence on lateral bending, the less pronounced oscillation of angular velocity, the intermittent nature of crawling, and the inter-individual variability of a number of locomotion-related parameters across different larvae.

We analyzed the trajectory of each single larva by tracking the forward velocity v of a posterior midline point (Figure 2A) and the head orientation angular velocity ω (Figure 2E). The former illustrates that forward locomotion indeed consists of consecutive steps (strides) that are characterized by the alternation of increased (green) and decreased (red) values of v. The latter confirms that lateral bending occurs both, during crawl-runs (weathervaning) and during stationary crawl-pauses (headcasts). The individual trajectory depicted in Figure 2 can be seen in full length in Figure 2—video 1, while the short epoch depicted in the inset is shown in Figure 2—video 2.

To characterize both types of oscillations we first perform a Fourier analysis of v and ω on each of 100 larva tracks (Figure 3A). Across larvae we observe a robust peristaltic rhythm with a crawling frequency f_C normally distributed around 1.4 Hz in line with earlier results (Heckscher et al., 2012; Mantziaris et al., 2020; Thane et al., 2023). We also confirm previous reports that lateral bending manifests a slower, more variable rhythm f_T around 0.4 Hz (Wystrach et al., 2016). Then we detect all strides performed by an individual animal and verify their stereotypical structure in terms of stride duration, resulting body displacement, and the phase-dependence of v (Figure 2D).

In order to quantify the phase-dependence of angular motion we interpolate each detected stride into a 64-bin oscillation cycle of period 2π. The progression of forward velocity from head to tail during the peristaltic cycle is illustrated in Figure 3C (top). To analyze the angular metrics during the peristaltic cycle we consider the absolute angular change ignoring the left or right direction of a turn. This analysis reveals a robust phase dependence of the angular velocity ω and acceleration , and of the bend angle θ_b (Appendix 2-Figure 1A). Specifically, the angular velocity ω exhibits a smooth unimodal dependency on the stride period peaking during the early stride phase and just before the head forward velocity reaches its maximum (Figure 3C). Nevertheless, ω during crawling is lower than during pauses as shown when we normalize it by dividing with the average computed during the pause epochs (Figure 3C, bottom). This implies a phasic attenuating interference of the peristaltic cycle on the lateral bending mechanism (Figure 1B). A plausible mechanistic explanation featuring bodily interference of crawling and bending is suggested in the Discussion.

Larvae intermittently pause crawling before re-assuming it. This results in sequences of concatenated strides (runs or stridechains) that alternate with brief crawl-pauses (Figure 2C). We analyzed the stridechain length with respect to the number of strides and the pause duration in the experimental dataset. To this end we pooled both measures across 100 larvae to obtain the empirical distributions. Testing power-law, exponential, log-normal and Levy distributions revealed the highest quality of fit for the log-normal distribution for both parameters (Figure 3B).

Simulation of behavioral experiments

In this section we simulate increasingly complex behavioral experiments. Starting from stimulus-free exploration we advance to chemotactic navigation and finally to adaptive odor preference experiments. We simulate individual larvae using the model as calibrated in Appendix 2 and all model parameters are shown in Appendix 2-Figure 1A. Larval populations are constructed by pooling individual virtual larvae that behave independently as they move through the spatial arena and odorscape (Niewalda et al., 2014).

Exploration

In the first virtual experiment, we consider free exploration (Figure 4) and evaluate our simulation results against an empirical data set (see Materials and Methods). To this end we simulate a population of 200 virtual larvae during three minutes while exploring a stimulus-free environment on a non-nutritious substrate in a Petri-dish, mimicking the experimental lab conditions (Figure 4). Due to the absence of sensory input and food, the proposed intermittent coupled-oscillator locomotory model at the basic layer is sufficient to autonomously generate exploratory movement.

Figure 4.

Larvae are initially placed at the center of the dish, each with a random body orientation. Over the course of the experiment, both the simulated and the real larval population disperses in space (Figure 4A, Figure 4—video 1).

The quantitative comparison shows a good agreement between simulated and empirical data with respect to the radial dispersal from the initial position (Figure 4B). Importantly, the variability across individual animals is well captured by the variability between individual model instances as shown for the number of strides, the crawling time, and the total distance travelled by each larva (Figure 4C).

Chemotaxis

Chemotaxis describes the process of exploiting an odor gradient in space to locate an attractive or avoid a repelling odor source. An olfactory sensor (olfactor) placed at the front end of the virtual body enables active sensing during body bending and allows detection of concentration changes that modulate turning behavior accordingly. Sensory-driven behavior is enabled via modulation of the frequency and amplitude of the lateral oscillator (turner) by sensory feedback, while the peristaltic oscillator is not affected, as proposed in (Wystrach et al., 2016).

To assess the chemotactic efficiency of our coupled-oscillator model we reconstruct the arena and odor landscape (odorscape) of two behavioral experiments described in (Gomez-Marin et al., 2011). In the first, larvae are placed on the left side of the arena facing to the right. An appetitive odor source is placed on the right side. The virtual larvae navigate up the odor gradient approaching the source (Figure 5A,C,E), reproducing the experimental observation in Fig.1 C in (Gomez-Marin et al., 2011). In the second, both the odor source and the virtual larvae are placed at the center of the arena. The larvae perform localized exploration, generating trajectories across and around the odor source. (Figure 5B,D,F), again replicating the observation in Fig.1 D in (Gomez-Marin et al., 2011). Two sample simulations can be seen in Figure 5—video 1.

Figure 5.

Odor preference test

We simulate the odor preference paradigm as described in the Maggot Learning Manual (Michels et al., 2017). Larvae are placed at the center of a dish containing two odor sources in opposite sides and left to freely explore. The odor concentrations are Gaussian-shaped and overlapping, resulting in an odorscape of appetitive and/or aversive opposing gradients. After 3 minutes the final situation is evaluated. The established population-level metric used is the olfactory preference index (PI), computed for the left odor as where N_l and N_r is the number of larvae on the left and right side of the dish while N is the total number of larvae.

The extend of olfactory modulation on the turning behavior is determined by the odor-specific gain G. As this is measured in arbitrary units, we first need to define a realistic value range that correlates with the behaviorally measured PI. We perform a parameter-space search independently varying the gain for left and right odors and measuring the resulting PI in simulations of 30 larvae. The results for a total of 25² gain combinations within a suitable range of G ∈ [−100, 100] are illustrated in Figure 6A. Simulation examples for one appetitive and one aversive odor are shown in Figure 6—video 1.

Figure 6.

In order to simulate larval group behavior in response to an associative learning paradigm we interface our behavioral simulation with the spiking mushroom body (MB) model introduced in (Jürgensen et al., 2024) (Figure 6C). It implements a biologically realistic neural network model of the olfactory pathway according to detailed anatomical data using leaky integrate-and-fire neurons (Jürgensen et al., 2021). The MB network undergoes associative plasticity at the synapses between the Kenyon cells and two MB output neurons as a result of concurrent stimulation with an odor and a reward signal. Both, odor and reward is simulated as spike train input to the receptor neurons and the reinforcement signalling dopaminergic neuron, respectively. The model employs two output neurons (MB₊, MB₋), representing a larger number of MB compartments associated with approach/avoidance learning respectively (Saumweber et al., 2018). The initially balanced firing rates between MB₊ and MB₋ are skewed after learning and encode the acquired odor valence (Owald et al., 2015; Owald and Waddell, 2015) here defined as:

We first trained the MB model via a classical conditioning experiment where, in each conditioning trial, it experiences an odor (conditioned stimulus, CS₊) in combination with a sugar reward during 5 min, following the standard training protocol in (Michels et al., 2017). 5 groups of 30 MB models undergo between 1 and 5 sequential conditioning trials (Weiglein et al., 2019). The resulting odor valence MB_out from each MB model was converted to an odor gain G via a simple linear transformation and used to generate a virtual larva (Figure 6B). Each population of 30 larvae was then tested in an odor preference simulation. The larvae were placed on a dish in presence of the previously rewarded odor (CS+) and a neutral odor in opposite sides of the dish (Figure 6—video 1), again following standard experimental procedures (Michels et al., 2017). To obtain robust results we replicated the experiment 100 times per population with a different random seed for a total of 600 simulations. The obtained preference indexes (PI) are illustrated in Figure 6D. The PI increase with increasing number of trials as well as its saturation resembles empirical observations (Weiglein et al., 2019). Note that the variability of the PI across the 100 simulations per condition is introduced solely by the behavioral simulation and resembles that seen across real experiments.

The current implementation only sequentially couples a trained MB model to be tested in a behavioral simulation. In the discussion we further elaborate on a possible extension featuring their closed-loop integration allowing for full behavioral simulations of both the training and the testing phase of the associative learning paradigm in a virtual environment.

Discussion

Experimental evidence for layered behavioral architecture

The proposed behavioral architecture is based on two underlying principles that justify our modular, hierarchically layered approach. Firstly, we attempt to horizontally structure it into behaviorally functional modules, each generating a well-defined behavior. Secondly, we vertically parse it into semi-independent layers, each under top-down modulation from the ones above but still capable of generating behavior independently of them. We summarize the neuroanatomical observations and behavioral experiments that support this approach.

The neural mechanisms underlying the three basic larval behaviors (Figure 1A, basic layer) have been extensively studied. Crawling is characterized by fairly stereotypical repetitive strides. Initial contraction of the head and tail segments driven by a ‘visceral piston’ mechanism is followed by a laterally symmetric peristaltic wave traversing neighboring segments longitudinally from back to front (Heckscher et al., 2012) (Figure 2B). Segmental central pattern generators (CPGs) coupled via intersegmental short- and long-range connectivity motifs involving also premotor neurons constitute the underlying neural circuitry (Kohsaka et al., 2019; Zarin et al., 2019; Mantziaris et al., 2020). Lateral bending results from asymmetric contraction of body musculature initiated at the thoracic segments (Berni, 2015). Feeding is generated via a network of mono- and multi-synaptic sensorimotor loops from enteric, pharyngeal and external sensory organs to motor neurons controlling mouth-hook movement, head-tilt and pharyngeal pumping (Miroschnikow et al., 2018b; Schoofs et al., 2024). A noteworthy facet of these behaviors is their autonomous generation with-out the need for descending control. This has been demonstrated by continued exploration even after brain ablation (Sims et al., 2019). Exploitation of a nutritious substrate requires all three basic behaviors as the larva constantly consumes food and re-positions its body. According to the layer-independence principle, we postulate that peripheral consummatory circuits at the sub-esophageal zone (SEZ) are also capable of autonomous exploitation.

Higher brain centers play a pivotal role in behavior modulation, with neurotransmitters like dopamine, serotonin, acetylcholine, and octopamine being key players in this intricate process (Berni et al., 2012; Zhang et al., 2013; Miroschnikow et al., 2018b; Malloy et al., 2019; Eschbach et al., 2020; Vogt et al., 2021). It has been suggested that the transition between exploration and exploitation (feeding) is acutely induced via dopaminergic signaling (Schleyer et al., 2020) while their long-term balance is regulated via hugin-mediated homeostatic neuromodulation (Schoofs et al., 2014). Identification of sensory pathways towards motor effector neuropiles further elucidates the role of interoception in behavioral modulation (Qian et al., 2018).

The neural mechanisms that underlie olfactory modulation of the basic locomotory behavior are also under intense investigation. Chemotactic approach and avoidance to innately valenced odors has been attributed to a predominantly innate pathway involving the antennal lobe (AL) and its direct projection to the lateral horn (LH), both in the larva (Schulze et al., 2015; Vogt et al., 2021) and in the adult fly (de Belle and Heisenberg, 1994; Strutz et al., 2014; Dolan et al., 2018). Modulation of behavior by learned odors strongly involves descending control by the MB in juvenile (Saumweber et al., 2018; Eschbach et al., 2020) and adult (Slater et al., 2015) stages. Both pathways have similar modulating effects on foraging behavior and are likely integrated in a premotor network downstream of the AL (Schleyer et al., 2015a; Eschbach et al., 2020). In the sensorimotor loop, descending pathways involving the LH control cessation of crawling, possibly triggering sharper re-orientation when navigating down-gradient, facilitating chemotaxis (Tastekin et al., 2018). Finally, the internal homeostatic state (e.g. starvation vs. satiation) is implicated in behavioral regulation via neurotransmitter release at multiple levels, including AL, LH and MB and neuropeptide expression (Vogt et al., 2021; de Tredern et al., 2024).

Previous computational models of larva locomotion

Larva locomotion has been modeled extensively. Generative models of peristaltic crawling feature neural and/or neuromuscular dynamics. Sequential neural activity patterns across the segmental body can be generated by longitudinally repeated CPGs of paired excitatory and inhibitory population rate units, possibly involving proprioceptive feedback (Gjorgjieva et al., 2013; Pehlevan et al., 2016). This relatively abstract CPG model has been elaborated into a connectome-based circuitry of premotor and motor neurons. Under tonic activation the model exhibits two functionally distinct, though structurally overlapping, recurring patterns of neuromuscular activity, responsible for forward and backward peristalsis respectively (Zarin et al., 2019). At the other end of the modeling spectrum the contribution of the visceral-piston mechanism to the peristaltic cycle has been assessed in a biomechanical model (Ross et al., 2015). A more elaborate neuromechanical model, based on segmental localized reflexes and substrate frictional forces and assuming empirically informed axial and transverse oscillatory frequencies, has been shown to generate forward and backward crawling without the need for any neural activation pattern (Loveless et al., 2019).

Lateral bending has also been captured in statistical (Davies et al., 2015) or generative models (Wystrach et al., 2016; Loveless et al., 2019; Loveless and Webb, 2018). In the context of free-exploration it has been shown to be a byproduct of chaotic body neuromechanics underlying peristalsis (Loveless et al., 2019). In the context of chemotaxis it has been attributed to a distinct oscillatory process, autonomous (Wystrach et al., 2016) or semi-autonomous to crawling (Davies et al., 2015). In the former case oscillation is driven by mutual inhibition between excitatory-inhibitory circuits. In the latter, bending behavior is dissected into low-amplitude weathervaning while crawling and high-amplitude headcasting during crawl-pauses, an approach essentially equivalent to an attenuation of the lateral oscillation amplitude due to crawling. Autonomous exploration and chemotaxis can be generated by an integration of the neuromechanical and the independent lateral oscillator models (Loveless and Webb, 2018).

Apart from the rare occasion where such models implement mutually exclusive mechanistic hypotheses and are indubitably incompatible to each other, in most cases they are indeed complementary, overlapping, nested or disconnected in terms of the generative mechanisms they aim to capture and could potentially coexist under a broader control system. In this context our unifying architecture for larval behavior in which any partial model can be positioned can be valuable for modelers and roboticists, interested in a behavior-based synthetic approach.

Intermittent coupled-oscillator model for realistic locomotion

Each of the above models could be adjusted so that at minimum it generates the 2D translational and angular motion of a virtual body and could therefore populate the basic locomotory layer of the behavioral architecture. In this study we propose such a model, assembled in a synthetic approach by distinct modules, which either extend previously proposed locomotory models or integrate novel findings derived from our analysis of kinematic parameters. Given a simplified bisegmental body, a simple oscillator under frequency-regulating tonic activation is adequate for generating recurrent strides, efficiently summarizing the complex underlying neural and neuromechanical dynamics (Heckscher et al., 2012; Mantziaris et al., 2020). Concerning angular motion, the previously introduced lateral oscillator model (Wystrach et al., 2016) meets the requirements and can therefore be coupled to the forward oscillator.

Our proposed model contributes two novel features. First, the intermittent nature of crawling as transitions between runs and pauses (Figure 3B). And, second, the peristaltic cycle-phase dependent attenuation of angular motion (Figure 3D). By combining these two features, the two bending behavioral modes termed weathervaning and head-casting are naturally generated via a phasic coupling between the two oscillators.

Behavioral intermittency

Larval locomotion is intermittent meaning that crawling runs are transiently intermitted by brief pauses. Transitions between these two behavioral states occur autonomously in stimulus-free conditions as during free exploration. Traditionally, in the context of movement ecology, intermittency has been studied in the spatial regime by characterizing the distributions of run distances and turn amplitudes occurring during brief stationary reorientation events (pauses). Power-law distributed runs in line with Levy-walk theoretical models (Günther et al., 2016; Sims et al., 2019) and diffusion-like kinematics have been reported (Klein et al., 2017). On the other hand, the turn amplitude distribution diverges from the original Levy-walk-predicted uniform distribution, as it is highly skewed towards small amplitudes, even if only significant reorientation events are taken into account (Sims et al., 2019). Moreover, the speed-curvature power-law relationship has been disputed (Zago et al., 2016; Marken and Shaffer, 2017; Zago et al., 2017; Marken and Shaffer, 2018). Regarding the temporal dynamics of intermittency, the duration distribution of these pauses is commonly neglected in traditional Levy-walk literature as they are usually characterised via the amplitude distribution of their resulting turning events. A recent analysis of larval tracks reported that the duration of pauses follows a power-law while that of activity bouts a log-normal distribution (Sakagiannis et al., 2020), partly in line with findings in adult-fly studies (Ueno et al., 2012; Reynolds et al., 2015). In our current study we found log-normal best fits for both pause duration and stridechain length. The disagreement over the pause duration might be attributed to the different timescale assessed in the two studies. Contrary to the high-resolution 180-sec tracks in this study, those analyzed previously lasted at minimum 1024 sec and up to 1 hour, allowing for the detection of longer pauses and imposing the necessity to fit over 4 orders of magnitude bypassing the apparent drop around 10 sec, also seen in Fig.3 of (Sakagiannis et al., 2020).

Concerning computational modeling, the contribution of behavioral intermittency to locomotion has not been adequately appreciated. Candidate generative models can either simply sample statistical distributions (Sims et al., 2019) or feature a generative mechanism that yields state transitions. Stochastic state-transitions have been included in a model of larva exploration yielding exponentially distributed epochs of both runs and stationary headcasts (Davies et al., 2015). At a mechanistic level, a recent study presented a simple binary-neuron model exhibiting state transitions between power-law and non power-law regimes via self-limiting neuronal avalanches and proposed a plausible underlying mechanism that explains initiation/cessation of crawling (Sakagiannis et al., 2020). All these attempts can be considered instantiations of a behavioral intermittency module (Figure 1B) controlling cessation and re-initiation of crawling, central for generating epoch transitions. In the present study, we chose stochastic sampling from the empirically fitted parametric distribution models.

Crawl-bend interference

Crawling includes mouth hook motion. Specifically, the first phase of a crawling stride consists of concurrent forward motion of head and tail segments, aided by a ‘visceral pistoning’ mechanism that generates forward displacement of the gut. Subsequently, the mouth hooks anchor the head to the substrate so that the second phase of peristaltic motion can drag all other segments forward as well, completing the stride (Heckscher et al., 2012). With respect to turning, it is still debated whether individual turns should be considered as discrete reorientation events that are temporally non-overlapping with crawling bouts (Sims et al., 2019), or whether lateral bending occurs in an oscillatory fashion generating turns both during crawling (weathervaning) and during pauses (headcasts) (Wystrach et al., 2016; Thane et al., 2019). The latter is supported by detailed eigenshape analysis confirming that larvae rarely crawl straight, rather forward locomotion is always accompanied by continuous small amplitude lateral bending (Szigeti et al., 2015). It follows that crawling does not exclude bending, rather the two strictly co-occur.

Crawling and bending partially recruit the same effector neural circuitry and body musculature at least at the level of the thorax. Peristaltic motion during crawling includes sequential symmetric bilateral contraction of all segments while bending occurs due to asymmetric unilateral contraction of the thoracic segments. This partial effector overlap could result in interference between the two processes. Indeed here we report a phase-dependent attenuation of angular velocity (Figure 3C). Attenuation reaches a minimum at a specific phase of the cycle, closely preceding the increase of head forward velocity (Figure 3C:inset). This coincides with the stride phase when the head stops being anchored to the substrate and is therefore free to move laterally. When applying a phase-dependent Gaussian attenuation kernel on angular velocity we managed to accurately reproduce the empirical relation (Appendix 2-Figure 1C).

A reasonable hypothesis would then be that the asymmetric thoracic contraction generating lateral bending becomes easier when the head is not anchored to the substrate therefore during a specific phase interval of the stride cycle. We propose that crawling phasically interferes with lateral bending because of these bodily constraints. A consequence of this hypothesis is that the amplitude of turns generated during crawl-pauses (headcasts) is larger in comparison to those generated during crawling (weathervaning) because during pauses the crawling interference to lateral bending is lifted. It is this phenomenon that dominates the description of larva exploration as a Levy-walk with non-overlapping straight runs and reorientation events, where weathervaning is neglected (Günther et al., 2016; Sims et al., 2019). Nevertheless, it has been included in a previous stochastic model of larva exploration, where it has been treated as entirely distinct to headcasts, occurring during crawl-pauses, via the application of differential constraints on both the angular velocity and the resulting turn amplitude (Davies et al., 2015). By implementing a coupled-oscillator locomotory model we avoid such dual treatment of headcasts and weathervaning.

Olfactory learning in closed loop behavioral simulations

We have reproduced the results of a basic associative learning experiment in the fruit fly larva (Schleyer et al., 2018; Jürgensen et al., 2024) by the open-loop simulation of classical conditioning trials and subsequent closed-loop behavioral simulation during the memory retention test for individual virtual larvae (Figure 6B-D). This modeling approach can be extended in multiple ways. First, the larva demonstrates a number of additional learning capabilities such as differential conditioning (Schleyer et al., 2011; Schleyer et al., 2015b; Schleyer et al., 2018), extinction learning (Felsenberg and Waddell, 2019; Lesar et al., 2021), and relief learning (Saumweber et al., 2018; Weiglein et al., 2019; Gerber et al., 2014; König et al., 2018). Interfacing neural network simulations with candidate circuit and synaptic mechanisms of plasticity with our behavioral model allows to directly compare virtual and empirical behavioral experiments, both at the level of individual and group assays. Second, while information about odor concentration is provided through active sensing, simultaneous input from a feeder module could activate the dopaminergic reward pathway required for synaptic plasticity in the mushroom body. At the same time, simulation of food uptake and energy expenditure will regulate the agent’s energy homeostasis. This would further allow realistic foraging scenarios with food depletion and competition among animals.

Closing the loop from active sensing to associative memory formation and behavioral control requires to synchronize a (spiking) neural network at the adaptive layer with the sensory module (reactive layer), and the locomotory and feeding modules at the basic layer (Figure 1A). This will enable the simulation of virtual larvae experiencing spatial and temporal dynamics in a virtual environment or on a robotic platform (Helgadottir et al., 2013), and it will allow to test model hypotheses on sensory-motor integration and to infer predictions for experimental interventions such as optogenetic stimulation (Saumweber et al., 2018) or genetic manipulation (Saumweber et al., 2011; Michels et al., 2011; Widmann et al., 2016; Springer and Nawrot, 2021).

Materials and Methods

We first describe the experimental dataset and the software package used in this study. Then we explicitly describe each of the computational modules that comprise the proposed behavioral architecture. The metrics used throughout the analyses of both real and simulated datasets are described in Appendix 1.

Dataset description

The larva-tracking dataset was obtained by M.Schleyer and J. Thoener at the Leipzig Institute of Neurobiology. It consists of 31 experimental groups. Each group of ∼ 30 third-instar larvae (Canton S) was placed on agarose-filled Petri dishes of 15 cm diameter with no particular stimuli and video-filmed from above at a framerate of 16 Hz for 3 minutes. During video-tracking, 12 points are detected along the longitudinal axis of each larva. For the present study groups were pooled together in a single population and timepoints of detected collisions have been excluded. A subset consisting of the 200 most complete, uninterrupted tracks was selected. The x-y coordinate timeseries have been filtered with a first-order butterworth low-pass filter with a cutoff frequency of 2 Hz in order to decrease tracking-related noise but retain the behaviorally relevant crawling frequency of f_C ≃ 1.4 Hz. The effect of inadequate and excessive filtering is illustrated in Video 3.

Software package and code availability

All data processing, data analyses and model simulations were performed using our freely available python package Larvaworld (https://pypi.org/project/larvaworld/), a behavioral analysis and simulation platform for Drosophila larva. In Larvaworld simulated and empirical data are treated indistinguishably, meaning that the exact same analysis pipeline and behavioral metrics are applied to both. Concerning modeling, the behavioral architecture we propose here comprises the backbone for the construction, extension, configuration and fitting of behavioral models in Larvaworld. Moreover, the intermittent coupled-oscillator model as introduced in the present manuscript is available for simulations alongside other preconfigured models. Extensive documentation can be found at https://larvaworld.readthedocs.io.

Inter-individual variability in virtual populations

The dominant methodological approach in animal behavioral modeling yields generative models that capture the average behavior of a group of animals. Each model therefore strives to generate behavior resembling that of an idealized average animal, neglecting inter-individual variability. To challenge this, we contrast it with a group-level modeling approach that preserves individuality within the population. Here, the generative model is instantiated by a group of non-identical animats, with parameter combinations drawn from a fitted multivariate normal distribution, preserving their pairwise correlations. We first measured a number of endpoint parameters across a population of 200 larvae and fitted a multivariate Gaussian distribution. A subset of three of these parameters, body-length l, crawling frequency f_C, and scaled stride displacement is shown in Figure 7. The univariate and bivariate empirical distributions are illustrated in blue, while the bivariate projections are shown in red.

Figure 7.

Module definition

The building blocks comprising the behavioral architecture are described as separate modules of defined input and output. The modular architecture only specifies the required placeholders and remains agnostic to the specific module implementations. Nevertheless, in what follows, the specific modular implementation for the proposed coupled-oscillator locomotory model will be described alongside the general module-placeholder definition. Calibration of each module on the empirical dataset and parameter specification is described in detail in Appendix 2. The parameters of the final model configuration are shown in Table 1.

Table 1.

Bisegmental body

The virtual body is the architecture’s interface with the simulated environment. It moves through space under the control of the motor effectors comprising the basic locomotory layer of the architecture, and thereby repositions the sensors introduced in the intermediate reactive layer. For the simplest case of 2D locomotion it is therefore imperative to define at least the forward velocity v and the angular orientation velocity ω. Any locomotory model dynamically generating these two parameters is adequate for an abstract point-body or a single-segment body implementation.

For the proposed coupled-oscillator model we choose a bisegmental body implementation, additionally featuring the bending angle (posture) θ_b between the front and rear segment. Following (Wystrach et al., 2016), the body is modeled as a torsional analog of the mass-spring damper model. Torque angularly accelerates an inertia I against angular damping that resists motion and viscoelastic forces that resist deformation (lateral bending). The original mass-spring damper model and its torsional analogue are defined by the equations:

We introduce simplified coefficients for angular damping z = ^Z/_I in sec⁻¹ and bend deformation resistance k = ^K/_I in sec⁻² as the dimensions of inertia ML² cancel out. Similarly, the external torque-per-unit-of-inertia is now in angular acceleration units This drive is generated by the turner module and will be described below. Additionally we introduce a crawl-phase dependent suppression of angular motion during crawling described by the coefficient c_CT (t) = c_CT (ϕ_C), which will also be described below.

In its original implementation (Wystrach et al., 2016) the torsional body model deliberately neglects two aspects of the real turning behavior of the larva. First, there is no distinction between bending ω_b and orientation ω angular velocities. Second, there is no correction of the bending angle θ_b due to forward motion. We tackle the first via the bisegmental body so that ω_b between the front and rear vector is distinct from the front vector’s ω. Regarding the second, we introduce a simple linear bending-angle correction as the rear vector is aligned to the front vector’s orientation during forward motion, according to the equation:

where d is the linear displacement during a timestep, and θ_b is the bending angle before and after correction, l is the body length and c_b = 1 is the bend correction coefficient.

Concerning the forward motion, the simple kinematic implementation of the crawler module directly generates forward velocity v without taking into account any biomechanical dynamics.

Crawler module

The crawler module generates forward velocity v, under a continuous (tonic) activation signal I_C and an additional intermittent initiation/cessation signal generated by the intermittency module. Velocity generation involves three parameters : the larval body-length l, the scaled displacement per stride and the crawling frequency f_C. The latter linearly reflects the tonic input I_C, which is kept constant during the short-duration simulations. For an individual virtual larva and f_C are set during initialization, where are the mean and standard deviation of .

In the present study crawling behavior is modeled as an oscillatory process. Each oscillation generates a cycle of forward velocity v increase-decrease resulting in displacement of the larva along the axis of its front-segment, simplistically modeling the result of exactly one peristaltic stride. After termination of a stride a new is sampled. An analytically tractable curve is fitted to the average empirical velocity curve measured during strides by setting all 5 crawler-related parameters to the median empirically measured value (Appendix 2-Figure 1B):

where ϕ_C is the instantaneous phase of the crawler oscillation iterating from 0 to 2π during an oscillatory cycle, is the maximum scaled velocity during the cycle and is the phase where the maximum velocity occurs.

Turner module

The turner module is defined as the generative model of a torque-like output A_T (t) under a continuous activation signal I_T (t). The (considered non-dimensional) output is scaled to angular acceleration by a coefficient c_T in sec⁻² and applied to the body as external bending drive (Equation 1):

which is then applied to the body yielding the instantaneous angular velocity ω(t) and angular acceleration bending the virtual body laterally. It has been previously suggested that an underlying oscillatory process drives alternating bending to the left and right side (Wystrach et al., 2016). We empirically confirmed this by detecting a slow,variable rhythm f_T ≈ 0.4 on the angular velocity timeseries of individual larvae. Based on this observation we adopt an oscillatory approach for the alternating lateral bending and extend the lateral oscillator model described in (Wystrach et al., 2016). It consists of two mutually inhibiting components each having an excitatory and inhibitory neuron (E_L&C_L vs E_R&C_R). The spike-response of a neuron to a cumulative input x is given by R(x, h) where h is the half-response threshold. The neuromodulatory activation I_T (t) excites both components but also affects the gain g(I_T) and time constant τ_H (I_T) of an additional adaptation component H_x of each neuron. The left and right components quickly settle in antiphase, while adaptation ensures that periodic transitions occur. Under the baseline activation the oscillations occur at a dominant frequency f_T . Perturbations of this external drive cause transient changes in both amplitude and frequency because both, up to transient loss of oscillation. This feature is exploited during olfactory modulation. The resulting oscillator activity is defined as the instantaneous difference in the firing rates A_T = (E_L − E_R).

Calibration of the turner module parameters in order to achieve a resulting angular motion that is in agreement to the empirical observations requires concurrent specification of the angularmotion related angular damping z, bend deformation resistance k and torque scaling coefficient c_T . The calibration process is described in detail in Appendix 2.

Crawler-turner coupling

We couple the crawler and turner oscillators by imposing a crawler-phase dependent attenuation c_CT (ϕ_C) on the angular velocity ω (Equation 1). During the entire stride-cycle ω is scaled by a baseline suppression coefficient . Our kinematic analysis revealed that this baseline attenuation is partially lifted during a phase interval of the stride-cycle. The maximum additional relief is defined by the relief coefficient . Concerning the mode of transition from to and back to , we study two implementations, each defined by an additional parameter. In the “square” mode there is an acute transition to maximum relief and back to baseline for a specific interval of the cycle . In the “phase” mode the transition is described as a Gaussian kernel reaching maximum relief at phase when crawl-induced angular attenuation is minimum and the anterior body is maximally free to bend laterally, as shown in the equation:

The parameter set is optimized to match the average empirically measured angular motion observed during the stride cycle. The optimization process used is described in Appendix 2. The optimal Gaussian kernel is shown in Figure 3.

Intermittency module

We define the intermittency module as a placeholder for any model capable of generating transitions between runs and pauses (or equivalently runs and headcasts) (Figure 1B). During pauses, the input to the crawler module I_C = 0 and therefore forward velocity v = 0 while during runs this is unaffected. For the specific study we implement a statistical model where the duration for both the run and pause states is sampled at each initiation from a distribution that has been fitted to the experimental data. In the case of pauses the duration is normally measured in time units. For runs we employ an equivalent metric and therefore measure them as number of consecutive crawling strides (stridechains). The fitted distributions are shown in Table 1.

Olfactory sensor

Olfaction is introduced in the intermediate layer of the behavioral architecture allowing chemotactic behavior. The olfactory sensor is located at the front end of the virtual larva therefore any reorientation and/or displacement influences sensory input. As in (Gomez-Marin et al., 2011) we assume that olfactory perception A_O relates to changes in odor concentration C according to the Weber-Fechner law, meaning that ΔA_O ∼ ln ΔC. We further add a decay term that slowly resets A_O back to 0. The rate of change is given in Equation 8 where c_O = 1 is the olfactory decay coefficient, G_i is the gain for odor i and C_i the respective odor concentration:

Concerning how the perceived olfactory stimulation A_O modulates exploratory behavior we adopt the mechanism proposed in (Wystrach et al., 2016) (Figure 5). According to this model, the turner activation I is perturbed from its baseline value within a suitable range :

Mushroom body module

The underlying spiking network model of the Drosophila larva olfactory pathway and mushroom body (MB) implements the model published by (Jürgensen et al., 2024). It consists of conductance-based leaky integrate-and-fire neurons. It encompasses 21 olfactory receptor neurons (ORNs) and the same number of projection neurons (PNs) and local interneurons (LNs) that receive input from the ORNs in a one-to-one manner (Figure 6C). LNs form inhibitory connections with the PNs. In the MB, each of the 72 Kenyon cells (KC) receive excitatory input from 2-6 random PNs. Feedback inhibition to the KCs is provided via the single GABAergic anterior paired lateral neuron (APL). All KC provide excitatory input to two MB output neurons (MBONs). The two MBONs represent two different types of output neurons that exist in the larval MB and either represent approach or avoidance behavior. Plasticity at the KC :: MBON synapses facilitates the association of odors with reward or punishment (associative learning). KC :: MBON synapses employ a two factor learning rule: Pre-synaptic activation by an odor stimulus (action potential of the KC) triggers an exponentially decaying eligibility trace e(t), which determines the window of opportunity for synaptic change. Additional neuromodulatory input r(t) from one of the two DANs will lead to a reduction in synaptic strength (Figure 6C) proportional to e(t) · r(t). The acquired imbalance between the outputs of the two MBONs (behavioral bias) represents the association of odors with rewards/punishments. An initially balanced interaction of excitatory and inhibitory feedback components from MBONs onto both DANs encodes the models’ learning history in form of a prediction error. The concept of prediction error suggests that associative learning about stimulus A is proportional to the difference between the reinforcement currently received with stimulus A minus the reinforcement predicted by stimulus A (previously learned association) (Rescorla, 1972). In vertebrates an implementation of this error signal has been demonstrated in dopaminergic neurons (Schultz, 2015, 2016) and a similar mechanism is proposed for learning in the insect MB (Riemensperger et al., 2005; Springer and Nawrot, 2021; Bennett et al., 2021; Jürgensen et al., 2024). The time-resolved behavioral bias is used to compute gains for odors, which bias motor output of the virtual larvae.

Videos

Here we provide the videos cited in the text. All videos are derived from real or simulated experiments as visualized via the Larvaworld software package, after some additional editing.

Acknowledgements

This project received funding from the Ministry of Culture and Science of the State of Northrhine Westphalia through the science network ‘iBehave’ (https://ibehave.nrw/) within the program “Netzwerke 2021”, and from the German Research Foundation through the Research Unit ‘Structure, Plasticity and Behavioral Function of the Drosophila mushroom body’ (DFG-FOR 2705, grant no. 365082554, https://www.unigoettingen.de/en/601524.html). PS received a stipendship within the Research Training Group ‘Neural Circuit Analysis’ (DFG-RTG 1960, grant no. 233886668). We would like to thank Michael Schleyer and Juliane Thoener for providing the tracking dataset and for valuable discussions on the data analyses. We thank Bertram Gerber for valuable discussion and Anna Morozova for substantial support with the generation of the video and figure material.

Additional files

Video 1. Locomotory model for Drosophila larva. The function of the locomotory model in Figure 1B is illustrated by gradually integrating its 4 modules (crawler, turner, oscillator-coupling, crawling-intermittency). In each of the 6 videos the implemented modules are shown in the inset.

Video 2. Free exploration simulation. A population of 25 real (left) or virtual (right) larvae is placed on a dish and left to freely explore. The body of the real larvae has been simplified into 2 segments as described in Video 4.

Video 3. Filter selection. The effect of inadequate or excessive filtering of the empirical larva recordings is illustrated. The left video shows the jittery original recording while the effect of lowpass filtering at cut-off frequencies of 4 Hz, 2 Hz and 0.5 Hz is shown on the rest. Selection of an intermediate 2 Hz cut-off frequency eliminates the unrealistic jitter while preserving the behaviorally relevant crawling frequency.

Video 4. Bisegmental larva-body simplification. The first video shows the original larva body as recorded by the tracker. 12 points are tracked along its longitudinal axis defining 11 segments while 22 points constitute the body contour. In the second video the body contour is dropped. In the third video an artificial rectangular contour is added for each body segment. In the last video the body-midline is segmented into 2 segments. The absolute head orientation angle θ is preserved while the single bending angle between the 2 segments is defined as .

Figure 2—video 1. The full-length trajectory (A, E) colored according to angular and forward velocity.

Figure 2—video 2. The short track epoch depicted in the insets (A, E) colored according to forward and angular velocity.

Figure 4—video 1. The temporal course of dispersal for real (left) and virtual (right) larvae shown in (A).

Figure 5—video 1. Time course of the two simulated chemotaxis experiments.

Figure 6—video 1. Two odor preference simulations, one with an appetitive and one with an aversive odor source placed on the left side of the dish. A non-valenced odor source is placed on the right side.

Appendix 1

Metric definition

Body length

The instantaneous body-length of an individual larva fluctuates during crawling due to sub-sequent stretching and contraction. Its histogram is well fitted by a Gaussian distribution (data not shown). Therefore individual larva length l is defined as the median of the midline length across time (total length of the line connecting all 12 midline points). All spatial parameters including displacement and velocity can be scaled to this body-length, converting spatial units m or mm to dimensionless body-length units. Scaled spatial metrics are denoted by an additional° over the metric symbol.

Segmentation and angular metrics

To specify the body segmentation providing the most suitable contact/rotation point for the definition of the bending ω_b and orientation ω angular velocities we analyse their relationship in a subset of 40 larvae. Tracking of 12 midline points allows computation of the absolute orientation of 11 body-segments and the respective 10 angles θ₁₋₁₀ between successive body segments (Figure 2F). We define θ as the head-segment absolute orientation in reference to the x-axis because this defines the movement orientation of the animal. We ask how ω results from the bending of the body as this is captured by the 10 angular velocities ω₁₋₁₀. The regression analysis depicted in Figure 1B shows as expected that ω depends primarily on the front angular velocities while this dependence decays as we move towards the rear segments, in line with previous studies (Lahiri et al., 2011). Timeshift analysis also shows that the front 3 angles change concurrently while angles further down the midline are increasingly lagging behind (data not shown). The correlation analysis depicted in Figure 1C shows that the sum of the front 5 angular velocities best correlates to ω. In other words the cumulative body bend of the front 5 segments best predicts head reorientation. Therefore we define the reorientation-relevant bending angle as (Figure 2F). The remaining 5 angles between the rear body-segments can safely be neglected as they do not contribute to reorientation. This analysis results in a segmentation of the body in a front and a rear segments of length ratio 5 : 6. The segmentation process is demonstrated in Video 4.

Forward velocity

To define forward velocity we need to choose which midline-point is most suitable to track and which velocity metric to use for defining the start and end of a stride. To this end we perform stride annotation of 3-minute tracks of a population of 20 larvae using each of 24 candidate instantaneous velocity metrics, namely the velocities of the 12 points, the component velocities of the rearest 11 points parallel to their front segment’s absolute orientation and finally the centroid velocity. To compare the candidate metrics we compute the spatial cv_s and temporal cv_t coefficient of variation of the annotated strides for each larva to assess how variant their time duration and displacement is. We finally compute the mean and across individuals. In Figure 1A the spatiotemporal stride variance is shown for each candidate metric. We choose the metric that provides the minimal spatial and temporal stride variance, assuming that strides of an individual larva are more or less stereotypical in both duration and displacement (Heckscher et al., 2012). Our study reveals that the centroid velocity is the most suitable metric for stride annotation. All spatial metrics are therefore computed via this point’s displacement.

Appendix 1—figure 1.

Track tortuosity

Track tortuosity is quantified by the straightness index (S.I), a metric previously used in larva-track analysis (Sims et al., 2019), computed by advancing a fixed time window along the track and calculating at each point the ratio of the straight line distance to the actual distance travelled. This index, which varies from 0 (no movement) to 1 (straight line movement), can capture very well the complexity of the movement at various scales (set by the window time frame) throughout the track. As the time window decreases the smoothing effect is also reduced revealing increasing track details. Different window sizes have been used ranging from 2 to 20 seconds in order to capture both large-scale spatial trajectories and small-scale local movements. Tortuosity was computed for each larva through time and revealed changes in movement as larva alternated between straight line relocation, changes of direction and different degrees of tortuosity. Hence, a change in the S.I. along a larva’s track captures the magnitude of the change in movement pattern from intensive, area-restricted searching movements (higher tortuosity) to extensive, straighter line movements (lower tortuosity), and vice versa, across a wide range of spatial scales.

Epoch annotation

Strides (S) are annotated using the scaled forward velocity First we apply fourier analysis to detect the dominant crawling frequency within a suitable range 1 ⩽ f_C ⩽ 2.5 (Figure 3A). From this we derive the reference stride duration Then epochs are annotated under a number of constraints (Figure 2C):

• Each stride is contained between two local minima.

• The local maxima contained in the epoch needs to exceed a threshold .

• The duration of the epoch t needs to range within . This allows individual strides to temporally vary without overlapping so that adjacent strides can be concatenated in stridechains.

After stride annotation the displacement due to each individual stride is computed for each larva and divided by the larva’s body-length (Figure 2A-C). The individual distributions are well fitted by Gaussians (data not shown). Therefore are defined as the average and standard deviation of the scaled displacement per stride for each larva.

Crawling runs (R) are defined as uninterrupted sequences of successive strides, also termed stridechains (Figure 2C). Stridechain length is the number of concatenated strides in a run and is a discrete metric equivalent to the continuous crawling run duration. If the locomotory model under evaluation does not generate v oscillations and therefore strides can not be detected, runs are annotated using the plain forward velocity v timeseries as in a previous study (Sakagiannis et al., 2020). We first define a suitable threshold v_thr by detecting the minima of the pooled v bimodal distribution. Then we define runs as epochs where constantly v ⩾ v_thr. Pauses are then defined as epochs containing no strides (or equivalently not overlapping with runs) and during which (or v ≤ v_thr).

Turn epochs (T) are contained between pairs of successive sign changes of orientation angular velocity ω. Annotated epochs of left (T_L) and right (T_R) turns yield the respective turn angle amplitudes and as the absolute total change of orientation angle θ (Figure 2G), which can then be pooled into the overall absolute turn-angle amplitudes θ_T . Notably by this definition turns include both headcasts (H) occurring during crawl-pauses and weathervaning (W) occurring during runs. Turns happening exclusively during runs and pauses yield the turn-angle amplitudes of weathervaning θ_W and headcasts θ_H respectively.

Appendix 2

Model calibration

Here we describe in detail the calibration procedure for a locomotory model of an average idealized individual based on a reference dataset of 150 larvae tracked while freely exploring a Petri-dish over 3 minutes. The pipeline consists of 3 initial steps, which can be performed independently of each other and set the configuration parameters of the crawler, intermittency and turner modules respectively. In a final step that builds upon the results of the initial three, the parameters of the crawl-bend phase-locked interference are set. The first two steps are directly defined by the kinematic analysis results. The third and the final steps involve an optimization process to reach the parameter set best fitting the empirical data. Notably, along with the intrinsic parameters of a module, a number of additional module-related parameters might need to be defined, such as the body-length for the crawler module and the angular-motion relevant physical parameters for the turner module.

Appendix 2—figure 1.

For the hereby described average-larva model no noise is introduced to any module input or output and no parameter variability is allowed across the population of the generated virtual larvae. The sole source of stochasticity then is due to the random distribution sampling processes operated during the behavioral simulations. For the analysis of the generated datasets the exact same pipeline is used as for the empirical dataset. The configuration of the final locomotory model is illustrated in Table 1.

Forward motion

For the Crawler module a set of 5 parameters need to be defined, along with a default body-length. This parameter-set is adequate to dynamically generate realistic forward velocity v oscillations as defined in Equation 3. All parameters are set to the median value across the empirical larva group. In Figure 1A the empirical distributions of the 5 crawler-module parameters are shown along with the median value selected for the average module.

Crawl intermittency

For the Intermittency module the sampling distributions for the run and pause epochs need to be defined. For both a continuous time-duration distribution is define while for the former an additional discrete distribution of number of concatenated strides per stridechain (run) is computed. In all cases, the pooled distribution of the given epoch across the entire dataset is approximated by power-law, exponential, Levy, log-normal and combined log-normal/power-law distributions. All distributions are truncated within an empirically observed range and are fully defined by their usual arguments. The best-fitting candidates from each class are sorted according to their Kolmogorov-Smirnov distance D_KS from the empirical distribution and the one exhibiting the minimum D_KS is selected for use in the model’s intermittency module (Table 1). A validation is shown in Figure 1D where the distributions of the generated epochs (blue) match the empirical ones (red).

Angular motion

Realistically calibrating the locomotory model in the angular domain is more demanding. The Turner module only generates a torque-equivalent output while the torsional-spring body model applies both angular damping and restorative force due to the body-bend angle (Equation 1). We need to calibrate the 3 involved physics parameters, namely the angular damping z in sec⁻¹, the bend deformation resistance k in sec⁻² and the torque coefficient c_T in sec⁻². Additionally the parameters of the turner module itself need to be specified. Here we study two implementations. The first is the above described neural oscillator, for which we need to define the baseline tonic input , the spike-response steepness coefficient n_T and the time constant τ_T . The second is a simple sinusoidal oscillator, for which we need to define the baseline amplitude and frequency f_T . In order to have an empirical reference for isolated angular-only motion while avoiding the effect of crawling on the angular motion and we select the detected pause epochs exceeding 3 seconds and derive the pooled distributions for 3 angular metrics, namely the bending angle θ_b, the angular velocity ω and the angular acceleration We run an optimization algorithm aiming to minimize the Kolmogorov-Smirnov distance between the simulated and empirical distributions for these 3 angular metrics. In each iteration a Turner module of a given configuration is simulated for 3 minutes while Equation 1 is applied to derive θ_b and ω. The obtained optimal turner and physics parameters are shown in Table 1.

Crawl-bend interference

For the crawl-phase dependent attenuation of angular motion we study two implementations, each having three parameters that define the attenuation coefficient c_CT (ϕ_C) during the stride-cycle (Equation 6). To achieve this we use a genetic algorithm (GA) optimization process to reach a best-fitting parameter-set for this module. Although optimization is only allowed to search the parameter space for these three parameters, each configuration is evaluated using a complete virtual larva model, combined with the optimal configurations of all other modules defined in the previous three calibration steps. During each iteration of the GA process a group of 30 model configurations is simulated for 3 minutes in free exploration conditions. The best 6 configurations are then mutated and/or combined in a novel generation. The evaluation function minimized is two-fold.Firstly the distributions of the three angular metrics derived from the simulated track, namely the bending angle θ_b, the angular velocity ω and the angular acceleration are compared to the empirically measured target distributions via the Kolmogorov-Smirnov distance as described above for the angular motion. Secondly, stride-cycle analysis is carried out for each simulated larva. Detected strides are interpolated in over n = 64 bins in a 2π cycle and direction-normalized by inverting the right-turning strides. The average stride-cycle curve is then computed for each of the three metrics and compared to the average empirically measured curve x via the evaluation metric defined below. The metric has been selected so that the different scales of angles, angular velocities and accelerations are normalized before being summed. The final evaluation metric to be minimized is defined as the cumulative error across the three angular metrics for both the distributions and the average stride-cycle curves.

Here we study two configurations for the crawl phase-dependent suppression of the angular motion, namely a smooth gaussian relief curve and an acute square transition as described in Equation 6.The target empirical and optimal model’s average curve for the angular velocity ω are shown in Figure 1B.

A sample simulation of the complete model is shown in Figure 1D. The Crawler-generated forward velocity oscillations (1st row) attenuate the angular motion in a phase-locked manner (2nd row) resulting in low angular velocity during runs and acute headcasts during pauses (3rd row). The bending angle is additionally restored to 0 during forward motion (4th row).