A neuromechanical model for Drosophila larval crawling based on physical measurements

Abstract

Background

Animal locomotion requires dynamic interactions between neural circuits, the body (typically muscles), and surrounding environments. While the neural circuitry of movement has been intensively studied, how these outputs are integrated with body mechanics (neuromechanics) is less clear, in part due to the lack of understanding of the biomechanical properties of animal bodies. Here, we propose an integrated neuromechanical model of movement based on physical measurements by taking Drosophila larvae as a model of soft-bodied animals.

Results

We first characterized the kinematics of forward crawling in Drosophila larvae at a segmental and whole-body level. We then characterized the biomechanical parameters of fly larvae, namely the contraction forces generated by neural activity, and passive elastic and viscosity of the larval body using a stress-relaxation test. We established a mathematical neuromechanical model based on the physical measurements described above, obtaining seven kinematic values characterizing crawling locomotion. By optimizing the parameters in the neural circuit, our neuromechanical model succeeded in quantitatively reproducing the kinematics of larval locomotion that were obtained experimentally. This model could reproduce the observation of optogenetic studies reported previously. The model predicted that peristaltic locomotion could be exhibited in a low-friction condition. Analysis of floating larvae provided results consistent with this prediction. Furthermore, the model predicted a significant contribution of intersegmental connections in the central nervous system, which contrasts with a previous study. This hypothesis allowed us to make a testable prediction for the variability in intersegmental connection in sister species of the genus Drosophila.

Conclusions

We generated a neurochemical model based on physical measurement to provide a new foundation to study locomotion in soft-bodied animals and soft robot engineering.

Background

Animal behaviour requires both neural circuits and the body. Specific neural circuits referred to as central pattern generators (CPGs) are capable of creating spatiotemporal activity patterns for behaviours [1,2,3,4,5]. CPGs regulate motor neuron activity and muscular contraction to generate force for movement. On the other hand, the body’s physical properties and constraints are also significant factors to guarantee animal movements [6,7,8,9]. Accordingly, the integration of neural circuits and body mechanics, referred to as neuromechanics, is one of the essential approaches to revealing motor control mechanisms [10,11,12,13,14,15].

Soft-bodied animals locomote by deforming their stretchable bodies [16]. Compared to animals with hard skeletons, soft-bodied animals possess high flexibility and can move adaptively in complicated environments [17]. This flexibility of soft-bodied animals suggests that output from neural circuits is not the sole determinant for locomotion. Animal’s physical properties, such as stiffness and viscosity, should also be involved in the dynamics of locomotion [9, 15] and have been measured experimentally [18,19,20,21,22,23]. However, it is still challenging to reproduce locomotion in soft-bodied animals by mathematical models based on physical measurements. To tackle this issue, we used Drosophila melanogaster (Drosophila, hereafter) larvae as a model and attempted to build a neuromechanical model describing their crawling behaviour.

Forward crawling is the most predominant mode in fly larval locomotion [24, 25]. By propagating segmental contraction from the posterior to anterior segments, larvae move forward [26,27,28], and neural circuits for crawling have been intensively examined [29,30,31,32,33,34,35,36,37]. Previous simulation studies have succeeded in building models describing the propagative nature of crawling behaviour qualitatively [10, 38,39,40]. In contrast, the physical properties of larvae remain to be studied, and quantitative reproduction of crawling kinematics has not been achieved yet.

Here, we present a neuromechanical model to describe fly larval crawling behaviour based on their locomotion kinematics and soft-bodied biomechanics. First, we recorded forward crawling in the third-instar larvae whose segmental boundaries were labelled by a fluorescent protein and extracted kinematics parameters of crawling. Then, using a tensile tester combined with optogenetics, we measured larval contraction force. Furthermore, larval viscoelasticity was measured by the tensile tester. The measurement indicated that the properties of the larval body were described better by the standard linear solid (SLS) model than the Kelvin-Voigt model previously used [10, 41]. By incorporating these physical parameters with a neural circuit model built previously [10], we established a neuromechanical model for larval locomotion. By optimizing the neural circuit model parameters, we succeeded in building a neuromechanical model that reproduced the kinematics parameter extracted from larva measurements. This model could also reproduce optogenetic studies reported previously. In addition, crawling in a low-friction condition was simulated, and its prediction was confirmed experimentally by analysing floating larvae. Perturbation analyses in our model predicted the importance of intersegmental connections in speed control, contrasting with previous studies [10, 40]. Based on our simulation results, we predicted the intersegmental connectivity in sister species in the genus Drosophila. Our model for larval crawling based on physical measurements provides a new approach to locomotion research in soft-bodied animals and soft robot engineering.

Results

Measurement of forward crawling behaviour in third-instar larvae

To reveal the physical mechanisms behind Drosophila larval crawling, we examined the kinematics of the behaviour and physical properties of the body. We used third-instar larvae because their size (3.53 ± 0.12 mm, n = 9 larvae) was large enough to measure physical properties as a soft material. While previous studies examined first- or second-instar Drosophila larvae and characterized the segmental kinematics [26, 27, 37, 42, 43], crawling behaviour in the third-instar larvae had not been investigated at a segmental scale but instead at a scale of the entire body [24, 25, 28, 32, 44]. Thus, we examined the segmental dynamics in freely moving third-instar larvae.

First, we recorded the position of the segment boundaries in freely crawling larvae. To reliably trace the segmental boundaries of larvae, we expressed a green fluorescent protein (GFP)-tagged coagulation protein Fondue, which accumulates at the muscle attachment sites [45]. We acquired time-lapse fluorescence images from the dorsal side of the larvae, where longitudinal muscles span single segments (Fig. 1A). From these fluorescence images, the segmental boundaries of the thoracic (T2–T3) and abdominal (A1–A8) segments were annotated (Fig. 1B; see the “Methods” section for details.) We examined the displacement of each segment boundary during forward crawling. As previously reported in the first- and second-instar larvae cases, sequential displacement of the segmental boundary from the posterior to anterior segments was observed in third-instar larvae (Fig. 1C). Four kinematic parameters were measured from the recording: stride length, stride duration, intersegmental delay, and speed (Fig. 1D). We obtained that stride length was 0.70 ± 0.05 mm, stride duration was 1.07 ± 0.12 sec, intersegmental delay was 0.09 ± 0.01 s, and speed was 0.64 ± 0.02 mm/s (n = 9 larvae).

[IMAGE OMITTED: SEE PDF]

Next, we quantified the dynamics of segmental length change. Plots of the segmental length showed the propagation of segmental contraction in crawling larvae (Fig. 1E, F). We measured three kinematic parameters (maximum segment length, minimum segment length, and contraction duration described in Fig. 1G) and obtained that maximum segmental length was 0.41 ± 0.03 mm, minimum segmental length was 0.21 ± 0.02 mm, and contraction duration was 0.44 ± 0.03 s (n = 9 larvae). To sum, we measured the seven kinematics quantities from segmental dynamics in freely crawling GFP-labelled third-instar larvae.

Consistent kinematic properties among segments in larval crawling

To establish a mathematical model for larval crawling, we examined whether we could assume that all the segments possess similar kinematic properties or not. To this aim, we compared the dynamics of segmental length change (defined in Fig. 1G) among all the segments (Fig. 2A–C). We noticed that the most posterior A8 segment had a shorter maximum length and smaller contraction range than the other segments (Fig. 2B). This observation reflected the fact that the A8 segment was a specialized structure and different from the others in the surface area of the body wall and the number of muscles [46]. In contrast, in the other segments, both the minimum and maximum lengths were comparable over the segments (Fig. 2B). Furthermore, contraction duration was also similar among the segments from T2 to A7 (Fig. 2C). Accordingly, most segments (T2–A7) exhibited similar kinematic dynamics in contraction during crawling, which allowed us to model the larva as a chain of segments with the same kinematic properties.

[IMAGE OMITTED: SEE PDF]

Having confirmed the similarity among the body segments, we inferred the effective number of segments. First, we adopted a segment model based on the “piston model” of larval peristalsis, where the motion of the head and tail were coupled [10, 26]. When a peristaltic wave reached the head, the tail was forced to move forward, which induced the contraction of the posterior segment A8 for the next stride of crawling. Then, we estimated the effective number of segments by assuming that the displacement of the virtual anterior end was the same as that of the tail. To evaluate the propagation velocity, we obtained intersegmental delay by calculating the cross-correlation of displacements between the most posterior segment boundary (y10) and the other segments (y0–y9) (Fig. 2D). The graph showed the propagation of segmental contraction from the posterior to the anterior segments. By plotting the intersegmental delay normalized by the stride duration, we noticed that the segmental contraction propagates at about a constant speed (Fig. 2E). This observation allowed us to model the larval crawling as the propagation of segmental contraction at a uniform speed. Based on this assumption, we inferred the effective number of segments of the whole body where the head and tail were coupled. To this aim, we fitted the intersegmental delay data by linear regression. The intersection points between the fitting lines and the horizontal line with a value of one on the vertical axis corresponded to the event when the propagation of segmental contraction reached the anterior end, and the subsequent propagation was initiated. From this calculation, we obtained 11.2 ± 0.2 as the effective number of segments (Fig. 2F, n = 9 larvae). Accordingly, the results indicated that the crawling larva could be modelled by a chain of 11 identical segments (Additional file 1: Fig. S1).

Measurement of larval contraction force

To establish a quantitative physical model for larval crawling, we attempted to measure the contraction force and soft material properties of larvae which enabled us to derive a Newtonian description of crawling by linking force to motion. We measured contraction force with the tensile tester (Fig. 3A). The tensile tester allowed us to measure and control the strain and stress of a sample material. We hooked an intact larva with insect pins and loaded it to the tensile tester. We measured the spontaneous contraction force in larvae and obtained values ranging from 1.4 to 2.7 mN (n = 11 spontaneous contraction events by three larvae, Fig. 3B, E). Although these values reflected the range of spontaneous contraction force, the larvae could not exhibit crawling locomotion during the measurement since their head and tail were hooked. Accordingly, we instead attempted to measure the maximum force that larvae could generate for crawling. To this aim, we adopted optogenetics to activate motor neurons reliably (Fig. 3C). We expressed a light-sensitive cation channel protein Channelrhodopsin2 [47] in all of the larval type-I motor neurons [48] and illuminated each larva with blue light (455 nm, 5.7 nW/mm2) to induce the contraction of all the body wall muscles. Since we needed to measure the muscular force evoked by the central nervous system (CNS), we expressed Channelrhodopsin2 in motor neurons instead of the muscles. The contraction force measurement with the optogenetic activation showed that the force ranged from 1.6 to 6.7 mN (n = 44 optogenetic stimulations applied to seven larvae, Fig. 3D, E). The variation in the measurement would be caused by the differences in some uncontrollable conditions, including hooking states at the head and tail. Consequently, we adopted the maximum measured force of 6.7 mN as the maximum force (\({F}_{M_{\mathrm{max}}}\)) used in our neuromechanical model.

[IMAGE OMITTED: SEE PDF]

Measurement of elastic and damping properties of larvae

Next, we analysed the passive properties of the larval body, elasticity, and viscosity, with the tensile tester (Fig. 4A). The stretchable body wall and hydrostatic skeleton of fruit fly larvae were typical components in soft-bodied organisms [49]. Previous studies reported material properties of soft-bodied animals, including the caterpillar Manduca sexta [18] and the earthworm Lumbricus terrestris [50, 51]. Although previous theoretical analyses of fly larval locomotion assumed that each segment was equivalent to a pair of a spring (a linear elastic component) and a damper (a linear viscous component) [10, 39, 41], the mechanical properties of Drosophila larvae were not measured experimentally. Accordingly, we measured the viscoelastic properties of fly larvae by methods used in material science.

[IMAGE OMITTED: SEE PDF]

First, we analysed the range of extension within which we could assume linearity between the extension and tensile force of the larval body (Fig. 4B). To analyse the relationship between the extension and tensile force of larvae, we stretched the larvae at a constant speed (Fig. 4C) and plotted tensile force against the extension (Fig. 4D). To test the linearity, we used linear regression. Depending on the range of extension, the goodness of fitting varies: linear regression fitted better in the range from 0 to 0.4 mm extension (Fig. 4E, left) than in a longer range from 0 to 1.0 mm (Fig. 4E, right). To find a suitable extension range with good linearity, we calculated the coefficient of determinants, which was an indicator of the predictability of the data by linear regression (Fig. 4F). The plot showed that extensions up to 0.4 mm or 0.6 mm exhibited high coefficients. Based on this observation, we concluded that we could assume linearity in the stretch ranging from 0 to 0.6 mm.

Next, we measured the viscoelasticity of larvae. Viscoelastic properties of soft materials could be acquired by measuring the relationship between stress (force loaded in the material) and strain (deformation of the material) [52]. To obtain the viscoelastic properties of larvae, we conducted the stress relaxation test, one of the standard dynamic mechanical tests [52]. In this test, a soft material sample was quickly extended to a certain length, and then the stress was being recorded over time (Fig. 4G). Based on the linearity test above, we extended the third-instar larvae by 0.4 mm. The result of the stress relaxation test showed an exponential decay that was approaching a non-zero plateau at a later time (Fig. 4H, left). To extract the physical properties from the stress relaxation data, a suitable physical model to describe the larval body was required. Among the general mechanical structure models of linear combinations of springs and dampers, we found that the Maxwell model and the Kelvin-Voigt model, which were used previously to model larval musculature [10, 39], could not fit the experimental curve well since either exponential decay or residual elastic force was missing in these models (Additional file 2: Fig. S2). In contrast, the standard linear solid (SLS) model, which combined the Maxwell model and a Hookean spring in parallel [52] (Fig. 4I), fitted the experimental results well. Thus, we fitted this curve with the general SLS model (Fig. 4H, left). Based on the SLS model, we obtained the two spring constants and one damping coefficient for the whole body of intact larvae by the extension of 0.4 mm (\({k}_1^{\mathrm{whole}}\) = 3.7 ± 0.7 mN/mm, \({k}_2^{\mathrm{whole}}\) = 5.3 ± 1.1 mN/mm, cwhole = (2.4 ± 0.5) ×102 mN s/mm, from n = 6 larvae, Fig. 4J–L). Consistent values were obtained by the extension of 0.6 mm (\({k}_1^{\mathrm{whole}}\) = 2.5 ± 0.4 mN/mm, \({k}_2^{\mathrm{whole}}\) = 6.1 ± 1.1 mN/mm, cwhole = (2.6 ± 0.2) ×102 mN s/mm, from n = 5 larvae, Fig. 4H right J–L), which further supported the linear property of larvae. The two spring constants did not have a strong correlation (r = 0.49, p = 0.13, n = 11; Pearson’s correlation coefficient) suggesting that they could not be reduced to a single spring constant. To sum, we obtained the viscoelastic parameters of larvae, two spring constants and one damping coefficient, by the stress relaxation test. It should be noted that in larval crawling, each segment exhibited contraction instead of extension which we induced in the measurement. However, we could not measure the viscoelasticity of larvae during contraction since larvae were bent instead of compressed when pushed from the both ends. Accordingly, we assumed larvae possess the same linear viscoelastic properties both in extension and compression. It should also be noted that in larval crawling, each segment contracted by about 0.2 mm (Fig. 2B) while in the stress-relaxation experiment, each segment was extended by about 0.04 mm (0.4 mm/11 segments). There could be some non-linearity in the viscoelastic properties of the larval body. However, in our model, we adopted the viscoelasticity measured with 0.4-mm extension as the approximate values for our mathematical model of larval crawling.

A neuromechanical model for larval crawling referring to the physical measurement

We established a mathematical neuromechanical model based on the physical measurements described above. On the other hand, we obtained seven kinematic values characterizing crawling locomotion (Fig. 5A) from the measurement of fly larvae (Figs. 1 and 2). Our aim was to build a model with a neural circuit and mechanical components to reproduce the observed results quantitatively. Regarding the mechanical part, based on the measurements, we estimated that the larva was a chain of eleven consistent segments (Fig. 2), each of which was described by the SLS model (Fig. 4 and Additional file 1: Fig. S1). These assumptions indicated that the viscoelastic parameters for single segments in the model (two spring constants k1 and k2 and one damping coefficient c) should be 11 times larger than the ones for the whole body (k1 = 40.7 mN/mm, k2 = 58.3 mN/mm, and c = 2640 mN s/mm). Furthermore, the force measurement indicated that the maximum muscular contraction force could be 6.7 mN (Fig. 3). We modelled segmental boundaries as masses that were pulled and pushed by the segmental muscles and dragged on the surface with friction (Fig. 5B). We measured the mass of larvae (1.14 mg, the 99% confidence interval is [1.12, 1.15], n = 70 larvae) and calculated the mass of one segment by dividing the mass of the whole body by the number of segments to obtain 0.1 mg/segment. The length of larvae was 3.5 ± 0.1 mm (n = 9). We used all these values in our mechanical model (green parameters in Fig. 5B). In contrast, for the neural network part, we adopted a neural circuit model built by Pehlevan et al. [10] based on Wilson-Cowan equations. This model was capable of generating propagation of waves from the posterior to anterior segments by interaction among excitatory and inhibitory neurons in the CNS, motor neurons, and sensory neurons [10]. Although the parameters in this model were suggested in the previous paper, we incorporated the neural circuit framework and tuned its parameters to build a neuromechanical model based on physical measurement to reproduce larval crawling.

[IMAGE OMITTED: SEE PDF]

Aiming to reproduce the measurement results on larval crawling, we adjusted 15 parameters in the neural circuit model (magenta parameters in Fig. 5B). The optimized set of the parameters (Additional file 3: Fig. S3) could reproduce the measurement results of the kinematics (Fig. 5C–F). To assess this observation quantitatively, we compared the seven kinematics parameters (Fig. 5A) between the measurement and simulation and found that all the parameters could be reproduced by the simulation to fit in the ranges of the measurement results (Fig. 5G–M). Accordingly, this neuromechanical model could quantitatively describe the segmental kinematics during forward larvae crawling.

The stress-relaxation results showed large variabilities on the order of 100% of the mean values (Fig. 4J–L). This variance could be caused by the uncertainty in the measurement or the variability in physical properties among larvae. Since we did not have independent measurement methods to obtain viscoelasticity, we could not determine which one was the case. However, when we perturbed the spring constants by the order of 100% of the optimized value, the crawling speed changed to be out of the measured range (Additional file 4: Fig. S4). This observation suggested that the viscoelasticity would be consistent between larvae, and the variabilities in the stress-relaxation test would be caused by uncertainty in the measurement.

We adopted μbf, the asymmetricity parameter between the forward and backward friction is 10, considering the fact that a larva had denticle belts on the ventral side of the body, providing asymmetric friction with the ground. To test the dependency of crawling on this asymmetricity, we perturbed μbf (Additional file 5: Fig. S5). The simulation results showed that crawling speed did not change even when μbf equalled 1, meaning that the friction was symmetric between forward and backward. This observation indicated that the asymmetricity in the friction along the body axis had a small role in generating forward crawling.

As for the maximum muscular contraction force, we adopted the maximum force observed by optogenetic stimulation. To test whether we could use smaller muscular forces to reproduce larval crawling, we conducted a perturbation analysis on the maximum muscular force in our neuromechanical model (Additional file 6: Fig. S6). The simulation results showed that crawling speed reduced as the muscular force decreased. Accordingly, this observation suggested that crawling speed was sensitive to the maximum muscular force, and the measurement of contraction force was crucial to building a neuromechanical model for larval crawling.

We compared parameters in our model based on physical measurements with those in the previous simulation by Pehlevan et al. [10]. In Pehlevan et al. [10], parameters were scaled by three values, the spring constant k, the equilibrium length of one segment L, and the relaxation time scale of the excitatory neurons τE. Accordingly, we substituted these parameters with our measurement results (k = 50 mN/mm, L = 0.35 mm, and τE = 35 ms) and obtained the absolute values of the parameters. By comparison, we found that the neural couplings and relaxation times in the CNS were almost consistent between our results and the assumption in Pehlevan et al. [10]. On the other hand, physical properties showed clear contrast between them. The damping coefficient in our result (c ~ 2640 mN s/mm) was about 400 times larger than the previous assumption [10] (c ~ 6 mN s/mm). When we adopted the damping coefficient 400 times smaller than our value, the crawling speed was affected to be doubled (crawling speed 0.69 mm/s in the optimized condition; 1.32 mm/s when the damping coefficient c was 400 times smaller than the optimized value of c, Additional file 4: Fig. S4C), and the maximum muscle contraction force in our results (FMmax = 6.7 mN) was 2.2 times smaller than the previous assumption [10] (FMmax = 15 mN). In addition, neural couplings between the CNS and either muscles or sensory neurons also exhibited differences: The coupling between the excitatory neurons and muscles in our model (wME = 540) was larger the previous assumption (wME = 1 in [10]), and the couplings between the sensory neurons and the CNS in our model (wES = wIS = 300) were also larger than the previous assumption (wES = wIS = 1.95 in [10]). To sum, our mechanical model had larger viscosity and smaller muscular force whereas the circuit model had larger couplings between the CNS and their targets in the peripheral structures than the assumption in the previous theoretical work.

The model reproduces previous observations in optogenetic experiments

To validate our model, we simulated two optogenetic perturbation studies reported previously and compared the results with the experimental observations [32, 53]. First, by optogenetic silencing of motor neurons in a few segments, propagation of contraction from the posterior to anterior segments was arrested, and all the segments were relaxed [53]. In the previous simulation study [10], although propagation was arrested by silencing the excitatory group in one segment, its neighbouring segment kept contracted instead of being relaxed. We silenced excitatory neurons in the A3 segment for 2 s in our neuromechanical model (Fig. 6A–F) and found that the silencing arrested the propagating wave (Fig. 6E). Furthermore, all the segments were relaxed by this optogenetic silencing (Fig. 6F), which is consistent with the experimental observation [53]. Inada et al. (2011) demonstrated that crawling was resumed after removing optogenetic silencing [53]. This phenomenon was not reproduced in our simulation: No wave was generated after optogenetic silencing (Fig. 6E, F), although waves were resumed after a short period of silencing (Additional file 7: Fig. S7). As pointed out in Pehlevan et al. [10], this would be due to the lack of such a mechanism in the neural circuit model. To sum, our model successfully reproduced the observation that the silencing excitatory neurons, particularly motor neurons, in a single segment arrested the propagation of peristaltic waves and relaxed all the segments.

[IMAGE OMITTED: SEE PDF]

Next, we tested the effect of activation of inhibitory neurons on crawling behaviour. A previous study reported that activation of inhibitory premotor interneurons PMSIs arrested wave propagation and caused the relaxation of the whole body [32]. Therefore, we activated inhibitory neurons in the A3 segment in the neuromechanical model (Fig. 6G). By this treatment, waves were arrested (Fig. 6H), and the whole segments returned to their equilibrium length (Fig. 6I). Accordingly, our neuromechanical model successfully reproduced the two optogenetics experiment results.

Measurement of forward crawling behaviour in a low-friction condition and validation of the model

Friction was a reaction force of a propelling force that drove larval movement. Although friction was necessary to move the centroid of larvae, it was unclear whether friction was required to generate the patterned segmental motion. In the previous theoretical study, two phenomena in segmental dynamics in a low-friction condition were simulated: Each segment contracted more, and intersegmental propagation delay became longer [10]. Our neuromechanical model tested the segmental dynamics in no friction conditions (Fig. 7A–H). The simulation results showed that the range of segment contraction did not change even in the no friction condition (Fig. 7C–F). Also, the duration of segmental contraction was not affected either (Fig. 7G). Furthermore, with respect to the wave propagation, the intersegmental delay was not changed in no friction condition (Fig. 7H). Accordingly, in our mathematical model, segmental dynamics was robust in the no friction condition.

[IMAGE OMITTED: SEE PDF]

To reveal the roles of friction on larval crawling, we measured the body segment kinematics in crawling larvae in a low-friction environment. We used a high concentration sugar solution (66% w/w sucrose) and floated third-instar larvae on it. To visualize the segment boundaries of larvae, we used tubP-Gal4, UAS-fondue::GFP fly larvae as described above (Fig. 1A). We compared the crawling motion between on the agarose substrate (Fig. 7I) and in the sugar solution (low friction environment) (Fig. 7J). Floating larvae could exhibit crawling behaviour (Fig. 7I), suggesting that the interaction with a substrate was not required to generate peristaltic motion. They could not move forward in the sugar solution, probably due to a lack of the propelling force that was the reaction force of friction or because of viscous resistance in the sugar solution. We also found that the number of waves was fewer in floating larvae than larvae on an agarose substrate (agarose 45.5 ± 3.5 crawls/min, n = 6 larvae; floating 6.8 ± 2.6 crawls/min, n = 6 larvae; p = 9.0 ×10−6, Welch’s t-test; Additional file 8: Fig. S8A). This observation implied that friction could be involved in the initiation of crawling. This phenomenon could not be observed in our neuromechanical model (Additional file 8: Fig. S8B) probably because our model did not include a mechanism to control the initiation of waves independent of the completion of waves. How each crawling wave was initiated and how the temporal gap between consecutive waves was regulated remained important open questions. Since our focus in this study was on the segmental and intersegmental dynamics, we analysed the kinematics of every single peristaltic wave.

First, we examined the segmental dynamics in the low friction condition and found that the range of segment contraction did not change in most of the segments in floating larvae (Fig. 7K, L). In addition, the contraction duration of each segment was not affected in the floating condition (Fig. 7M). These observations were consistent with our model results (Fig. 7E–G). In the posterior segments A6 and A7, however, there was a difference between larvae on agarose and floating larvae. The range of contraction at the posterior segments was reduced in floating larvae (Fig. 7L). Since forward crawling started at the posterior segments, this observation implied that friction would be crucial to initiate crawling. This implication was consistent with the observation of the reduction in crawling numbers in floating larvae (Additional file 8: Fig. S8A). Second, we analysed the intersegmental delay in the floating larvae. In most segments, the intersegmental delay was consistent between larvae on agarose and floating larvae (Fig. 7N). However, the intersegmental delay between T3 and A1 was increased in floating larvae. Since this phenomenon was observed in the limited anterior segment, it would be possible that other behaviour, including head-sweeping that occurs in the anterior segments, was affected in the low friction condition. To sum, the kinematics in larval crawling was almost similar between on the agarose substrate and in low-friction condition. This observation was consistent with the result of our neuromechanical model.

Contribution of sensory feedback and CNS

The CNS and proprioception were both involved in locomotion. The significant contribution of sensory feedback in larval crawling was shown experimentally [44]. This observation raised a hypothesis that proprioception could cause the propagation signal to the neighbouring segment through intersegmental sensory feedback [44]. A simulation based on this assumption showed that the model circuit with intersegmental feedback could normally generate propagation waves without an intersegmental connection in the CNS [10]. On the other hand, several key interneurons regulating larval locomotion had been identified [32, 54,55,56], suggesting that interneurons in the CNS should be involved in wave propagation. Aiming to reproduce the experimental observation of the involvement of the CNS in crawling, we tested the contribution of intersegmental sensory feedback and the intersegmental connection in the CNS to peristaltic motion by our neuromechanical model (Fig. 8). When the sensory feedback was silenced (wES = 0), the crawling speed was reduced (Fig. 8B, D, and H). This was consistent with the previous simulation result [10] and the experimental result [44]. When the intersegmental central connection was blocked (wEn = 0), the speed of crawling was decreased too (Fig. 8B, F, and I). This result was also consistent with observations that interneurons in the CNS were involved in crawling speed. While the speed was reduced in both cases, the neural mechanisms would be different. When sensory feedback was blocked, the activity of inhibitory neurons was attenuated (Fig. 8C, E). Meanwhile, when the intersegmental central connection was blocked, the activity of inhibitory neurons was enhanced (Fig. 8C, G). This observation suggested that sensory feedback and the CNS were involved in the speed control by distinct mechanisms. To examine their roles in detail, we perturbed wES and wEn from the optimal values. The speed of crawling was robust to the change of wES (Fig. 8H), whereas increase in wEn led to a drastic change in crawling speed (Fig. 8I). In particular, crawling speed depended on wEn nonlinearly with a peak. To sum, this simulation result was consistent with the experimental observations that interneurons in the CNS should be involved in the regulation of crawling behaviour.

[IMAGE OMITTED: SEE PDF]

Next, we analysed the interaction between intersegmental sensory feedback and intersegmental connection in the CNS. We perturbed wEn and wES simultaneously and plotted the displacement of larvae in a two-dimensional diagram (Fig. 9A). The plot showed that crawling could be realized in the wide area of this parameter space. When wEn was larger than the value in the optimized condition (0.6), wEn played a dominant role in speed control compared to wES (Fig. 9B). On the other hand, in the regime of smaller wEn, both wEn and wES were involved in crawling speed. This result indicated the strength of central connectivity in the CNS should play a pivotal role in controlling the speed of crawling behaviour.

[IMAGE OMITTED: SEE PDF]

Contribution of excitatory and inhibitory neurons in crawling

Based on the observation suggesting a significant contribution of the CNS to crawling behaviour, we next examined the involvement of excitatory and inhibitory neurons in larval locomotion. To this aim, we perturbed intrasegmental connection weight: wEE for excitatory connections and wEI for inhibitory connections (Fig. 10A–G). Crawling speed decreased when blocked excitatory connections (Fig. 10B, D). Since the traces of excitatory and inhibitory neurons were different from those in blocking intersegmental excitatory connections (Figs. 8E and 10E), the speed reduction mechanism should be different between silencing the intrasegmental excitatory connections and intersegmental excitatory connections. On the other hand, when the activity of inhibitory neurons was suppressed, larvae could not exhibit crawling (Fig. 10F). In this case, all the neurons exhibited tonic hyperactivity due to a lack of inhibition (Fig. 10G). To further reveal the involvement of wEE and wEI, we perturbed them in the range where larvae could exhibit crawling. The plots showed that speed was related to wEE and wEI in non-monotonical manners (Fig. 10H, I). This observation indicated that both excitatory and inhibitory neurons controlled crawling speed.

[IMAGE OMITTED: SEE PDF]

Next, we analysed the interaction between excitatory and inhibitory neurons. Displacement of larvae with different values of wEE and wEI was plotted two-dimensionally. The diagram indicated that when wEE was large and wEI was close to zero (the upper right side of the diagram), larvae could not generate crawling (Fig. 11A). On the other hand, when wEE was small (weak excitatory connection) or wEI was small (strong inhibitory connection) that corresponded to the lower-left region of the diagram, and crawling speed was reduced. Interestingly, as long as wEE and wEI were balanced (the diagonal region of the diagram), crawling speed was consistent. The deviation in speed when wEE and wEI were perturbed was within the range of variability observed experimentally (Fig. 11B). These results indicated two points: The balance between excitatory and inhibitory connections within segments was critical to generating peristaltic motion, and, as long as they were balanced, crawling speed would not change by perturbation in intrasegmental connection weights.

[IMAGE OMITTED: SEE PDF]

Prediction in the neural network in sister species of Drosophila melanogaster

A comparison between the contribution of intersegmental (Fig. 9B) and intrasegmental (Fig. 11B) connections to crawling speed suggested that intersegmental connections in the CNS had dominant roles in regulating crawling speed. We noticed that this observation could allow us to predict neural circuits for fly larval crawling. A recent study reported an inter-specific variation in crawling speed among eleven species in the genus Drosophila [57]. This study showed that some sister species crawled faster than Drosophila melanogaster, and the others locomoted slower. We attempted to make a prediction about the intersegmental connection in the CNS of these species based on their crawling speed. To predict the property of neural circuits in these sister species, we overlaid the speed data of the species on the graph showing speed dependency on wEn (Fig. 12). To compensate the difference in experimental conditions between this study and Matsuo et al. [57], the speed of each species was normalized by the speed of Drosophila melanogaster. To analyse the dependency on wEn, the dependency of wES was averaged out from Fig. 9B. The result showed that all of the eleven horizontal lines demonstrating the crawling speed of sister species crossed to the speed-wEn relation curve (Fig. 12). This plot suggested that some species showing faster crawling, Drosophila pseudoobscura for example, should have strong intersegmental connections in the CNS while those crawling slower such as Drosophila mojavensis, should have weaker intersegmental central connections. This prediction would be testable by interspecific comparison based on the anatomy of their central nervous system or physiological recording.

[IMAGE OMITTED: SEE PDF]

Discussion

In this work, we built a neuromechanical model based on the physical properties of fly larvae. The model successfully described crawling behaviour quantitatively. First, we quantified segmental dynamics in larval crawling and obtained seven parameters that characterize the kinematics of peristalsis (Figs. 1 and 2). Then, using a tensile tester, viscoelasticity and muscular force of larvae were measured (Figs. 3 and 4). These results suggest that the larval body could be modelled as a chain of mass and SLS modules (Fig. 4), distinct from the previous modelling for larvae [10, 39, 41]. By incorporating the material properties and optimizing the parameters in the neural circuit, our model succeeded in reproducing larval crawling quantitatively (Fig. 5). Furthermore, this model could describe the observation in optogenetic studies (Fig. 6) and crawling in a low-friction condition tested in this study (Fig. 7). Perturbation analyses indicated the importance of intersegmental connections in the CNS, contrasting with previous studies (Figs. 8 and 9). Based on this observation, we predicted the intersegmental connections in the CNS in sister species in the genus Drosophila (Fig. 12).

One striking finding from our model contrasts with previous studies: the significant contribution of intersegmental connections in the CNS to crawling speed. Although several interneurons involved in speed control have been reported [32, 54,55,56], previous simulation studies were not consistent with these observations. Our model had larger viscosity and smaller contraction forces (Additional file 3: Fig. S3) than previous parameter regimes, implying that each segment is hard to move. The large weight of neuromuscular connection (wME) and that of sensory feedback (wES) compensated for this difficulty (Additional file 3: Fig. S3). With the strong interaction between the CNS and peripheral organs (muscles and proprioceptors), the CNS could play a significant role in controlling crawling speed (Fig. 8I). This finding demonstrates the significance of intersegmental connection in the CNS. It also highlights the importance of mechanical properties in the larval body to understand the dynamics of larval locomotion.

We measured the viscoelastic properties of larvae by the stress-relaxation test. It should be noted that there were some limitations in this measurement. First, we obtained these passive properties of the body by measuring unanaesthetized larvae. We analysed the data without any spontaneous contraction in the stress-relaxation test. However, there remained a possibility that infinitesimal muscular tension was induced by proprioceptive feedback upon stretching. By blocking action potentials in motor neurons while keeping spontaneous vesicle release, it would be able to suppress feedback from proprioception to measure purer viscoelasticity. Second, we approximated the viscoelasticity for the contraction in the range of 0.2 mm, which larvae showed in crawling, by the viscoelasticity for the extension in the range of 0.04 mm, which we adopted in the stress-relaxation test. In real larval crawling, asymmetricity between extension and contraction and nonlinearity in the viscosity would affect the kinematics. By incorporating these factors, it would be able to reproduce the kinematics more quantitatively.

The model in this study could describe forward crawling, the most frequent behaviour in larvae. Recent neuroscience studies have revealed numerous neural circuit modules involved in distinct aspects of behaviour, including speed control [32], bilateral coordination [35], intrasegmental coordination [36], intersegmental coordination [31], backward crawling [58, 59], turning [60], escaping [61], and sensory input-guided navigation [34, 62]. By integrating these circuit modules, it would be possible to establish a neuromechanical model that reproduces multiple and natural larval behaviours. Furthermore, the neuromechanical model in this study can potentially serve biomimetics. In recent decades, more and more soft robots, endowed with new capabilities relative to the traditional hard ones, have been designed to exhibit complex movements [63,64,65,66]. Taking Drosophila larvae as a prototype, bionic structures can be established with high dexterity to explore the unstructured environments. By harnessing perturbation analyses on the neuromechanical model, the locomotion ability of the soft larval robots could be tuned and optimized for required applications.

Conclusions

We established a neuromechanical model based on physical measurements of fly larvae to describe larval crawling motion. The model succeeded in quantitatively reproducing the segmental dynamics of larval locomotion measured experimentally in this work. The model also reproduced previous optogenetic observations. In addition, this model predicted a regular peristaltic motion in a low-friction condition, and this property was confirmed experimentally in this study. Furthermore, the model predicted an essential contribution of intersegmental neural connections in the central nervous system, contrasting with a previous report. This hypothesis provided a testable prediction for the diversity in intersegmental connectivity in sister species in the genus Drosophila. This study provided a new foundation to investigate crawling kinematics in soft-bodied animals and soft robot engineering.

Methods

Fly strains

We used the third-instar larvae in all the experiments. To quantify segment dynamics, we used tubP-Gal4, UAS-fondue::GFP (Bloomington #5138, #43646), which labels the segment boundaries [45]. In the viscoelasticity measurements and the spontaneous contraction force measurement, we used a wild-type strain Canton S. For the optogenetic experiment, OK6-Gal4 [67] and UAS-ChR2[T159C]::YFP [31] were used.

Image acquisition and segment boundary annotation

Third-instar larvae of tubP-Gal4, UAS-fondue::GFP were gently washed in deionized water to remove residual food from the body surface. Then, individual larvae were placed on a flat agarose stage (1.5% agarose). Behavioural videos were captured by a charge-coupled device (CCD) camera (XCD60, Sony, Japan) mounted on an Olympus stereomicroscope (SZX16, Olympus, Japan). Images were acquired at 30 Hz. We used the Fiji software [68] to manually annotate the right and left ends of every segment boundary (from the anterior boundary of the T2 segment to the posterior boundary of the A8 segment; see Fig. 1A, n = 9 larvae, biological replicates). The midpoint of the right and left ends of a boundary was used as the position of the boundary.

Spontaneous contraction force and viscoelasticity measurement

Third-instar larvae in the feeding stage were washed with deionized water and dried with paper. Insect pins were inserted into the head and tail. The pin in the head was bent to form a loop to hook to a paper clip that was hung on the hook of the tensile sensing machine, while the pin in the tail was used to fix the body on the PDMS (Polydimethylsiloxane, Sylgard 184, Toray, Japan) silicone block, held by the tong on the SHIMADZU EZ-S platform with a 5-N load cell. The experimental equipment is shown in Fig. 3A (for the measurement of spontaneous contraction force) and Fig. 4A (for the stress relaxation test). The larval body was kept on the vertical axis for measurement. The baseline of force was calibrated by values measured before applying external elongational force. All the experimental procedures were performed at room temperature.

In the measurement of the spontaneous contraction force, we recorded eleven contraction events from three larvae that are biological replicates (Fig. 3E).

During the stress relaxation test to measure viscoelasticity, we applied a constant strain of 0.4 mm, about 10% of the body length (3.53 ± 0.12 mm, n = 9 larvae, biological replicates), to the larvae (n = 6 larvae, biological replicates). Then, the stress decreased until it reached a plateau after a while, as shown in Fig. 4H. To fit the stress relaxation curve to mechanical analogues, we adopted the Maxwell, Kelvin-Voigt, and standard linear solid (SLS) models. These models give the relationship between joint force F and displacement ∆L. Detailed functions are described as follows:

$$\frac{F}{c}+\frac{1}{k}\frac{dF}{dt}=\frac{d\Delta L}{dt}\kern0.75em \mathrm{Maxwell}\ \mathrm{model}$$

$$F=k\Delta L+c\frac{d\Delta L}{dt}\kern0.75em \mathrm{Kelvin}-\mathrm{Voigt}\ \mathrm{model}$$

$$F+\frac{c}{k_2}\frac{dF}{dt}={k}_1\left(\Delta L+c\frac{k_1+{k}_2}{k_1{k}_2}\frac{d\Delta L}{dt}\right)\kern0.75em \mathrm{SLS}\ \mathrm{model}$$

where k (or k1and k2 in the SLS model) and c are the elastic and damping constants, respectively, and t is time. The total displacement of 0.4 mm was realized by pulling the larval body at 1 mm/min. During the stress relaxation experiments, this elongation time (24 s) is much shorter than that for relaxation (576 s), and thus, the displacement is regarded as the step function L(t) = L0H(t − t0) and initial force is F(0) = 0N, where L0 = 0.4 mm, and H(t) is the Heaviside step function. In this case, the corresponding stress relaxation functions in these models are described as follows:

$$F=k{e}^{-\frac{k}{c}\left(t-{t}_0\right)}{L}_0H\left(t-{t}_0\right)\kern0.5em \mathrm{Maxwell}\ \mathrm{model}$$

$$F=k{L}_0H\left(t-{t}_0\right)+c{L}_0\delta \left(t-{t}_0\right)\kern0.5em \mathrm{Kelvin}-\mathrm{Voigt}\ \mathrm{model}$$

$$F(t)=\left[{k}_1+{k}_2{e}^{-\frac{k_2\left(t-{t}_0\right)}{c}}\right]{L}_0H\left(t-{t}_0\right)\kern0.5em \mathrm{SLS}\ \mathrm{model}$$

where t0 ≥ 0. By fitting these curves to the stress-relaxation measurement data, we obtained the spring constants and damping coefficients. We used Python 3.7 for the curve fitting.

Contraction force measurement with optogenetics

For optogenetic activation of motor neurons, we used the OK6-GAL4, UAS-ChR2 line [32]. The early third-instar larvae were selected and put into ATR (all-trans-retinal) containing yeast paste, the concentration of which was 1 mM. These larvae were reared at 25 °C in the dark for 1 day [69]. Afterwards, the third-instar larvae were prepared to measure tension force as described above (the “Spontaneous contraction force and viscoelasticity measurement” section). Blue light-emitting diode (LED) light (455 nm, 5.7 nW/mm2, M455L3, ThorLab) was used to stimulate Channelrhodopsin2 (ChR2) expressed in motor neurons, which leads to the contraction of the larval body. In each stimulation, the blue light was applied for 2 s, and the evoked force was monitored. Seven larvae (biological replicates) were used in the measurement, and each measurement took two to 5 min. The optogenetically induced forces were measured by the differences between forces during illumination and no illumination (Fig. 3).

Larval crawling in a low-friction environment

We used a high concentration sugar solution (66% w/w sucrose) and floated third-instar larvae in it. The genotype we used was tubP-Gal4, UAS-fondue::GFP, which allowed us to mark the segmental boundary from T2 to A8 (Fig. 1A, B). The locomotion (n = 6 larvae, biological replicates) was recorded via an SZX16 fluorescent microscope (Olympus, Japan) with a × 1.25 object lens at 30 frames/s, and its trajectory was measured by Fiji [68].

Modelling

Body-substrate mechanics and modelling

The larval crawling stride consists of the piston and wave phases [26]. The piston phase constrains the larval body length to be almost constant, while the wave phase generates the propagation wave repetitively. To make it simple, we modelled the whole body of larvae as a chain of eleven segments.

We assumed that the segments in the model larva are entirely repetitive based on the observation in Figs. 1 and 2. Then, the body was modelled as a chain of the SLS units in series (Additional file 1: Fig. S1). The muscle groups were modelled as the tension actuator to accept efferent control from the neural circuit (Fig. 5B). \({F}_{M_i}\) is the muscular tension force in the ith segment to counteract the effect of body viscoelasticity and friction. The mechanics are described based on Newton’s second law as follows.

$$\left\{\begin{array}{c}m{\ddot{y}}_{-1}=m{\ddot{y}}_{10}={F}_{10}-{F}_0+{F}_{M_{10}}-{F}_{M_0}-{f}_{-1}\kern7em \\ {}m{\ddot{y}}_i={F}_i-{F}_{i+1}+{F}_{M_i}-{F}_{M_{i+1}}-{f}_i\kern1.5em i=0,\dots, n-2\kern3.25em \\ {}{y}_{-1}={y}_{10}+ nL\kern21em \end{array}\right.$$

$${F}_i+\frac{c}{k_2}\frac{d{F}_i}{dt}={k}_1\left({y}_{i-1}-{y}_i-L+c\frac{k_1+{k}_2}{k_1{k}_2}\left(\frac{d{y}_{i-1}}{dt}-\frac{d{y}_i}{dt}\right)\right)$$

$${\tau}_M{\dot{F}}_{M_i}=-{F}_{M_i}+{F}_{M_{\mathrm{max}}}{\sigma}_M\left({w}_{\mathrm{ME}}{E}_i-{\theta}_M\right)\kern2em i=0,\dots, n$$

$${\sigma}_M(v)=\frac{1}{1+{e}^{-{k}_Mv}}$$

where yi is the segmental boundary position (Fig. 1B and Additional file 1: Fig. S1), m is the mass, Fi is the viscoelastic force, \({F}_{M_i}\) is the muscular tension force, fi is the friction force (see below), τM is the time relaxation constant for the muscular tension, wME is the connection weight from the excitatory unit to the muscular actuator, θM is the threshold for the muscular tension, kM is the gain in function σM(v), and n is the total number of segments, equal to 11. y−1 is the position of the head, and y10 is the position of the tail. To model the piston phase, where the head and tail move concurrently, the posterior and anterior ends share the same velocity during crawling, as \({\dot{y}}_{-1}={\dot{y}}_{10}\) which is inferred from the third equation above. Each mass is subject to a joint force of the force from SLS modules, the muscular tension, and the friction force. The force from the SLS unit is affected by the segmental displacement, and the muscular tension force is regulated by the sigmoid function with gain kM.

As for the friction force, it exists during the interaction of the mechanical body with the substrate. We modelled that the joint force on a mass was zero when the speed of the mass equals zero, and the sum of the force from the SLS module and the muscular force was weak by giving the counterforce by the friction. In other cases, The friction was Ff in forward movement and μbfFf in backward movement. The ratio of forward to backward friction (μbf) was introduced to the model as directional asymmetric, considering the anterior-posterior polarity in denticle bands. Since the backward friction is larger than the forward friction, the ratio is more than one. Accordingly, the friction on the mass was represented as:

$${f}_i=\left\{\begin{array}{c}{F}_{i,\operatorname{ext}}\kern12em {\dot{y}}_i=0\ \mathrm{and}-{\mu}_{bf}{F}_f\le {F}_{i,\operatorname{ext}}\le {F}_f\\ {}{F}_f\left({\sigma}_f\left({\dot{y}}_i\right)-{\mu}_{bf}{\sigma}_f\left(-{\dot{y}}_i\right)\right)\kern2.5em \mathrm{otherwise}\kern9.25em \end{array}\right.\kern1.25em i=-1,0,\dots, n-1$$

$${\sigma}_f(v)=\frac{1}{1+{e}^{-{k}_fv}}$$

where fi, Ff, μbf (μbf > 1) individually represent the friction, forward friction, and ratio for forward-backward friction, respectively, and Fi, ext refers to the joint force of the force from the SLS unit and the muscular tension. As described above, when the mass is still and Fi, ext does not exceed the range of forward-backward friction, the friction force fi should be equal to Fi, ext to cancel out the total force (the joint force) on the mass. Otherwise, the friction is either Ff or μbfFf depending on the movement direction.

Neuromuscular dynamics and modelling

The framework for the neural circuit is depicted in Fig. 5B. The neural circuit model is based on the model in Pehlevan et al. [10]. Neural dynamics is realized by the activities of excitatory and inhibitory populations of neurons in each segment under the Wilson-Cowan model [70, 71]. In Drosophila larvae, most proprioceptive neurons were active when the segment was contracted [43]. Referring to the model in Pehlevan et al. [10], the feedback from the sensory receptor in our model provides signals for the inhibitory neuron within the same segment and the excitatory neuron in the anterior segment. In this case, it can promote both local relaxation and forward propagation waves. This connectivity is consistent with the “mission-accomplished” model [44].

Under these assumptions, the neural dynamics of neural circuits were described as follows:

$$\left\{\begin{array}{c}{\tau}_E{\dot{E}}_i=-{E}_i+{\sigma}_E\left({w}_{\mathrm{EI}}{I}_i+{w}_{\mathrm{EE}}{E}_i+{w}_{\mathrm{En}}{E}_{i+1}+{w}_{\mathrm{ES}}{S}_{i+1}-{\theta}_E\right)\kern1.75em i=0,\dots, n-2\kern0.75em \\ {}{\tau}_I{\dot{I}}_i=-{I}_i+{\sigma}_I\left({w}_{\mathrm{IE}}{E}_i+{w}_{\mathrm{IS}}{S}_i-{\theta}_I\right)\kern13.5em i=0,\dots, n-2\kern0.75em \\ {}{\tau}_E{\dot{E}}_{10}=-{E}_n+{\sigma}_E\left({w}_{\mathrm{EI}}{I}_{10}+{w}_{\mathrm{EE}}{E}_{10}+{w}_{\mathrm{En}}{E}_0+{w}_{\mathrm{ES}}{S}_0+\mathrm{ini}-{\theta}_E\right)\kern6.5em \\ {}{\tau}_I\dot{I_{10}}=-{I}_{10}+{\sigma}_I\left({w}_{\mathrm{IE}}{E}_{10}+{w}_{\mathrm{IS}}{S}_{10}-{\theta}_I\right)\kern18.75em \end{array}\right.$$

$${S}_i={\sigma}_S\left(L{\theta}_S-{y}_{i-1}+{y}_i\right)\kern13em i=0,\dots, n-1$$

$${\sigma}_E(v)=\frac{1}{1+{e}^{-{k}_Ev}},\kern2em {\sigma}_I(v)=\frac{1}{1+{e}^{-{k}_Iv}},\kern2em {\sigma}_S(v)=\frac{1}{1+{e}^{-{k}_Sv}}$$

where Ei and Ii are the mean firing rates of the excitatory and inhibitory population in the ith segment, respectively; τE and τI are the relaxation time constant of the excitatory and inhibitory units, respectively; wab is their synaptic connection weight from population b to population a; θE and θI are the activation threshold of the excitatory and inhibitory units, respectively; kE, kI, and kS are the sigmoid gain of excitatory, inhibitory, and sensory units, respectively; and Si is the feedback strength from the sensory receptor in the ith segment. The circuit includes excitatory connections via wIE and wME and inhibitory connections via wEI. When the segmental length becomes smaller than the threshold LθS, the sensory unit starts to be activated. To trigger the initial crawling, we introduced an external stimulus ini, a rectangular pulse (for 5 ms), on the excitatory population in the posterior terminal.

The values of the parameters are listed in Additional file 3: Fig. S3. All the simulation work is performed using stiff solver ode15s in MATLAB R2021a.

Availability of data and materials

All data generated or analysed during this study are included in this article and its supplementary information files. The segment kinematics data generated in this study are available from the Figshare repository with the DOI of 10.6084/m9.figshare.19289594 [72].

Abbreviations

CPGs:

Central pattern generators

GFP:

Green fluorescent protein

T:

Thoracic

A:

Abdominal

CNS:

Central nervous system

SLS:

Standard linear solid

CCD:

Charge-coupled device

PDMS:

Polydimethylsiloxane

ATR:

All-trans-retinal

LED:

Light-emitting diode

ChR2:

Channelrhodopsin2

References

Grillner S, El Manira A. Current principles of motor control, with special reference to vertebrate locomotion. Physiol Rev. 2020;100(1):271–320.

Ijspeert AJ. Central pattern generators for locomotion control in animals and robots: a review. Neural Netw. 2008;21(4):642–53.

Marder E, Calabrese RL. Principles of rhythmic motor pattern generation. Physiol Rev. 1996;76(3):687–717.

Marder E, Bucher D. Central pattern generators and the control of rhythmic movements. Curr Biol. 2001;11(23):R986–96.

Minassian K, Hofstoetter US, Dzeladini F, Guertin PA, Ijspeert A. The human central pattern generator for locomotion: does it exist and contribute to walking? Neuroscientist. 2017;23(6):649–63.

Gomez-Marin A, Ghazanfar AA. The life of behavior. Neuron. 2019;104(1):25–36.

Scheffer LK, Meinertzhagen IA. A connectome is not enough – what is still needed to understand the brain of Drosophila? J Exp Biol. 2021;224(21):1–9.

Miller LA, Goldman DI, Hedrick TL, Tytell ED, Wang ZJ, Yen J, et al. Using computational and mechanical models to study animal locomotion. Integr Comp Biol. 2012;52(5):553–75.

Tytell ED, Holmes P, Cohen AH. Spikes alone do not behavior make: why neuroscience needs biomechanics. Curr Opin Neurobiol. 2011;21(5):816–22.

Pehlevan C, Paoletti P, Mahadevan L. Integrative neuromechanics of crawling in D. melanogaster larvae. Elife. 2016;5:1–23.

Boyle JH, Berri S, Cohen N. Gait modulation in C. elegans: an integrated neuromechanical model. Front Comput Neurosci. 2012;6:1–15.

Ekeberg Ö, Grillner S. Simulations of neuromuscular control in lamprey swimming. Philos Trans R Soc B Biol Sci. 1999;354(1385):895–902.

Di RA, Stanev D, Armand S, Ijspeert A. Sensory modulation of gait characteristics in human locomotion: a neuromusculoskeletal modeling study. PLoS Comput Biol. 2021;17(5):1–33.

Ijspeert AJ, Crespi A, Cabelguen JM. Simulation and robotics studies of salamander locomotion. Neuroinformatics. 2005;3(3):171–95.

Nishikawa K, Biewener AA, Aerts P, Ahn AN, Chiel HJ, Daley MA, et al. Neuromechanics: an integrative approach for understanding motor control. Integr Comp Biol. 2007;47(1):16–54.

Elder HY. Peristaltic mechanisms. In: Aspects of animal movement [Society for Experimental Biology Seminar Series, vol5]. Cambridge: Cambridge UP; 1980. p. 71–92.

Alexander RM. Principles of animal locomotion. Princeton, New Jersey: Princeton University Press; 2006. p. 384.

Lin HT, Dorfmann AL, Trimmer BA. Soft-cuticle biomechanics: a constitutive model of anisotropy for caterpillar integument. J Theor Biol. 2009;256(3):447–57.

Dorfmann AL, Woods WA, Trimmer BA. Muscle performance in a soft-bodied terrestrial crawler: constitutive modelling of strain-rate dependency. J R Soc Interface. 2008;5(20):349–62.

Paterson BA, Marko Anikin I, Krans JL. Hysteresis in the production of force by larval Dipteran muscle. J Exp Biol. 2010;213(14):2483–93.

Backholm M, Ryu WS, Dalnoki-Veress K. Viscoelastic properties of the nematode Caenorhabditis elegans, a self-similar, shear-thinning worm. Proc Natl Acad Sci U S A. 2013;110(12):4528–33.

Gilpin W, Uppaluri S, Brangwynne CP. Worms under pressure: bulk mechanical properties of C. elegans are independent of the cuticle. Biophys J. 2015;108(8):1887–98.

Rabets Y, Backholm M, Dalnoki-Veress K, Ryu WS. Direct measurements of drag forces in C. elegans crawling locomotion. Biophys J. 2014;107(8):1980–7.

Berrigan D, Pepin DJ. How maggots move: allometry and kinematics of crawling in larval Diptera. J Insect Physiol. 1995;41(4):329–37.

Green CH, Burnet B, Connolly KJ. Organization and patterns of inter- and intraspecific variation in the behaviour of Drosophila larvae. Anim Behav. 1983;31(1):282–91.

Heckscher ES, Lockery SR, Doe CQ. Characterization of Drosophila larval crawling at the level of organism, segment, and somatic body wall musculature. J Neurosci. 2012;32(36):12460–71.

Lahiri S, Shen K, Klein M, Tang A, Kane E, Gershow M, et al. Two alternating motor programs drive navigation in Drosophila larva. PLoS One. 2011;6(8):e23180.

Fox LE, Soll DR, Wu CF. Coordination and modulation of locomotion pattern generators in Drosophila larvae: effects of altered biogenic amine levels by the tyramine β hydroxlyase mutation. J Neurosci. 2006;26(5):1486–98.

Hunter I, Coulson B, Zarin AA, Baines RA. The Drosophila larval locomotor circuit provides a model to understand neural circuit development and function. Front Neural Circuits. 2021;15:1–12.

Gowda SBM, Salim S, Mohammad F. Anatomy and neural pathways modulating distinct locomotor behaviors in Drosophila larva. Biology (Basel). 2021;10(2):1–30.

Kohsaka H, Zwart MF, Fushiki A, Fetter RD, Truman JW, Cardona A, et al. Regulation of forward and backward locomotion through intersegmental feedback circuits in Drosophila larvae. Nat Commun. 2019;10(1):2654.

Kohsaka H, Takasu E, Morimoto T, Nose A. A group of segmental premotor interneurons regulates the speed of axial locomotion in Drosophila larvae. Curr Biol. 2014;24(22):2632–42.

Clark MQ, Zarin AA, Carreira-Rosario A, Doe CQ. Neural circuits driving larval locomotion in Drosophila. Neural Dev. 2018;13(1):1–10.

Kohsaka H, Guertin PA, Nose A. Neural circuits underlying fly larval locomotion. Curr Pharm Des. 2017;23(12):1722–33.

Heckscher ES, Zarin AA, Faumont S, Clark MQ, Manning L, Fushiki A, et al. Even-Skipped+ interneurons are core components of a sensorimotor circuit that maintains left-right symmetric muscle contraction amplitude. Neuron. 2015;88(2):314–29.

Zwart MF, Pulver SR, Truman JW, Fushiki A, Fetter RD, Cardona A, et al. Selective inhibition mediates the sequential recruitment of motor pools. Neuron. 2016;91(3):615–28.

Zarin AA, Mark B, Cardona A, Litwin-Kumar A, Doe CQ. A multilayer circuit architecture for the generation of distinct locomotor behaviors in Drosophila. Elife. 2019;8:1–34.

Gjorgjieva J, Berni J, Evers JF, Eglen SJ. Neural circuits for peristaltic wave propagation in crawling Drosophila larvae: analysis and modeling. Front Comput Neurosci. 2013;7:1–19.

Ross D, Lagogiannis K, Webb B. A model of larval biomechanics reveals exploitable passive properties for efficient locomotion. Lect Notes Comput Sci (including Subser Lect Notes Artif Intell Lect Notes Bioinformatics). 2015;9222:1–12.

Loveless J, Lagogiannis K, Webb B. Modelling the neuromechanics of exploration and taxis in larval Drosophila. PLoS Comput Biol. 2019;15(7):1–33.

Paoletti P, Mahadevan L. A proprioceptive neuromechanical theory of crawling. Proc R Soc B Biol Sci. 2014;281(1790):20141092. https://doi.org/10.1098/rspb.2014.1092.

Berni J. Genetic dissection of a regionally differentiated network for exploratory behavior in drosophila larvae. Curr Biol. 2015;25(10):1319–26.

Vaadia RD, Li W, Voleti V, Singhania A, Hillman EMC, Grueber WB. Characterization of proprioceptive system dynamics in behaving Drosophila larvae using high-speed volumetric microscopy. Curr Biol. 2019;29(6):935–44.

Hughes CL, Thomas JB. A sensory feedback circuit coordinates muscle activity in Drosophila. Mol Cell Neurosci. 2007;35(2):383–96.

Green N, Odell N, Zych M, Clark C, Wang ZH, Biersmith B, et al. A common suite of coagulation proteins function in drosophila muscle attachment. Genetics. 2016;204(3):1075–87.

Bate M. The development of Drosophila melanogaster. Long Island, New York: Cold Spring Harbor Laboratory Press; 1993.

Pulver SR, Pashkovski SL, Hornstein NJ, Garrity PA, Griffith LC. Temporal dynamics of neuronal activation by channelrhodopsin-2 and TRPA1 determine behavioral output in Drosophila larvae. J Neurophysiol. 2009;101(6):3075–88.

Sanyal S. Genomic mapping and expression patterns of C380, OK6 and D42 enhancer trap lines in the larval nervous system of Drosophila. Gene Expr Patterns. 2009;9(5):371–80.

Trueman ER (Edwin R. The locomotion of soft-bodied animals. London: Edward Arnold; 1975. 200.

Quillin K. Ontogenetic scaling of hydrostatic skeletons: geometric, static stress and dynamic stress scaling of the earthworm Lumbricus terrestris. J Exp Biol. 1998;201(12):1871–83.

Quillin KJ. Kinematic scaling of locomotion by hydrostatic animals: ontogeny of peristaltic crawling by the earthworm Lumbricus terrestris. J Exp Biol. 1999;202(6):661–74.

Banks HT, Hu S, Kenz ZR. A brief review of elasticity and viscoelasticity for solids. Adv Appl Math Mech. 2011;3(1):1–51.

Inada K, Kohsaka H, Takasu E, Matsunaga T, Nose A. Optical dissection of neural circuits responsible for Drosophila larval locomotion with halorhodopsin. PLoS One. 2011;6(12):e29019.

Fushiki A, Zwart MF, Kohsaka H, Fetter RD, Cardona A, Nose A. A circuit mechanism for the propagation of waves of muscle contraction in Drosophila. Elife. 2016;5:1–23.

Yoshikawa S, Long H, Thomas JB. A subset of interneurons required for Drosophila larval locomotion. Mol Cell Neurosci. 2016;70:22–9.

Babski H, Surel C, Yoshikawa S, Valmier J, Thomas JB, Carroll P, et al. A GABAergic Maf-expressing interneuron subset regulates the speed of locomotion in Drosophila. Nat Commun. 2019;10(4796):1–17.

Matsuo Y, Nose A, Kohsaka H. Interspecies variation of larval locomotion kinematics in the genus Drosophila and its relation to habitat temperature. BMC Biol. 2021;19(1):1–21.

Hiramoto A, Jonaitis J, Niki S, Kohsaka H, Fetter RD, Cardona A, et al. Regulation of coordinated muscular relaxation in Drosophila larvae by a pattern-regulating intersegmental circuit. Nat Commun. 2021;12(1):1–14.

Carreira-Rosario A, Zarin AA, Clark MQ, Manning L, Fetter RD, Cardona A, et al. MDN brain descending neurons coordinately activate backward and inhibit forward locomotion. Elife. 2018;7:1–28.

Tastekin I, Riedl J, Schilling-Kurz V, Gomez-Marin A, Truman JW, Louis M. Role of the subesophageal zone in sensorimotor control of orientation in drosophila larva. Curr Biol. 2015;25(11):1448–60.

Ohyama T, Schneider-Mizell CM, Fetter RD, Aleman JV, Franconville R, Rivera-Alba M, et al. A multilevel multimodal circuit enhances action selection in Drosophila. Nature. 2015;520(7549):633–9.

Louis M. Mini-brain computations converting dynamic olfactory inputs into orientation behavior. Curr Opin Neurobiol. 2020;64:1–9.

Kim S, Laschi C, Trimmer B. Soft robotics: a bioinspired evolution in robotics. Trends Biotechnol. 2013;31(5):287–94.

Aguilar J, Zhang T, Qian F, Kingsbury M, McInroe B, Mazouchova N, et al. A review on locomotion robophysics: the study of movement at the intersection of robotics, soft matter and dynamical systems. Rep Prog Phys. 2016;79(11):110001.

Trivedi D, Rahn CD, Kier WM, Walker ID. Soft robotics: biological inspiration, state of the art, and future research. Appl Bionics Biomech. 2008;5(3):99–117.

Corucci F, Cheney N, Giorgio-Serchi F, Bongard J, Laschi C. Evolving soft locomotion in aquatic and terrestrial environments: effects of material properties and environmental transitions. Soft Robot. 2018;5(4):475–93.

Aberle H, Haghighi AP, Fetter RD, McCabe BD, Magalhães TR, Goodman CS. Wishful thinking encodes a BMP type II receptor that regulates synaptic growth in Drosophila. Neuron. 2002;33(4):545–58.

Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9:676–82.

Matsunaga T, Fushiki A, Nose A, Kohsaka H. Optogenetic perturbation of neural activity with laser illumination in semi-intact drosophila larvae in motion. J Vis Exp. 2013;77:1–5.

Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J. 1972;12(1):1–24.

Negahbani E, Steyn-Ross DA, Steyn-Ross ML, Wilson MT, Sleigh JW. Noise-induced precursors of state transitions in the stochastic Wilson–Cowan model. J Math Neurosci. 2015;5(1):1–27.

Sun X, Kohsaka H. Tracing of fly larval segments. figshare; 2022. https://doi.org/10.6084/m9.figshare.19289594.