Development of equation of motion deciphering locomotion including omega turns of Caenorhabditis elegans

eLife assessment

This useful study introduces a simple mechanical model of C. elegans locomotion that captures aspects of the worm's behavioral repertoire beyond forward crawling. While the kinetic model (ElegansBot) provides a compromise and starting point to help understand the mechanical components of C. elegans behavior, the claim that this work improves on extant mechanical models is incomplete, including modeling a 3-dimensional turning behavior with a 2-dimensional model without sufficient justification. In addition, the results of the application of the model to previously unstudied behaviors are primarily qualitative and do not produce new predictions.

https://doi.org/10.7554/eLife.92562.3.sa0

Significance of the findings:

Useful: Findings that have focused importance and scope

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Incomplete: Main claims are only partially supported

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Appendix 1
Data availability
References
Article and author information

Abstract

Locomotion is a fundamental behavior of Caenorhabditis elegans (C. elegans). Previous works on kinetic simulations of animals helped researchers understand the physical mechanisms of locomotion and the muscle-controlling principles of neuronal circuits as an actuator part. It has yet to be understood how C. elegans utilizes the frictional forces caused by the tension of its muscles to perform sequenced locomotive behaviors. Here, we present a two-dimensional rigid body chain model for the locomotion of C. elegans by developing Newtonian equations of motion for each body segment of C. elegans. Having accounted for friction-coefficients of the surrounding environment, elastic constants of C. elegans, and its kymogram from experiments, our kinetic model (ElegansBot) reproduced various locomotion of C. elegans such as, but not limited to, forward-backward-(omega turn)-forward locomotion constituting escaping behavior and delta-turn navigation. Additionally, ElegansBot precisely quantified the forces acting on each body segment of C. elegans to allow investigation of the force distribution. This model will facilitate our understanding of the detailed mechanism of various locomotive behaviors at any given friction-coefficients of the surrounding environment. Furthermore, as the model ensures the performance of realistic behavior, it can be used to research actuator-controller interaction between muscles and neuronal circuits.

Introduction

With only a few hundred neurons, Caenorhabditis Elegans (C. elegans) perform various behaviors such as locomotion, sleeping, reproduction, and hunting (Hall and Altun, 2008). The connectome structure among 302 neurons and 165 somatic cells of C. elegans was discovered by pioneering works (Cook et al., 2019; White et al., 1986). C. elegans is a cost-efficient and widely used model in neuronal research. Its small body size and minimal nutritional requirements contribute to its cost efficiency. The organism matures in a shorter period, about three days, compared to other model animals such as fruit flies or mice. Its transparent body allows for easy microscopic observation of its internal structures or artificially expressed green fluorescent proteins. Moreover, due to the hermaphroditic nature of C. elegans, offspring mostly share the same genotype as the parent, which simplifies the multiplication of the worm population for research purposes (Hall and Altun, 2008).

C. elegans bends its body with a sinusoidal wave pattern when moving forward or backward. The driving force for this movement comes from the difference between perpendicular and parallel frictional forces, which it experienced from a surrounding environment. This thrust force pushes the worm along the ground surface with which the worm contacts (Berri et al., 2009; Boyle et al., 2012; Hu et al., 2009; Niebur and Erdös, 1991). Even if a worm has a sinusoidal modulation generated inside it, it has difficulties in forward and backward locomotion if it does not feel the difference in frictional forces from its surroundings.

Mechanical simulators of rod-shaped animals such as C. elegans (Boyle et al., 2012; Niebur and Erdös, 1991), fish (Ekeberg, 1993), and snakes Hu et al., 2009 have been used in various studies. These simulators demonstrate how the activities of muscle cells are represented as behavioral phenotypes, which are determined by signals from a neuronal circuit simulator (Boyle et al., 2012; Ekeberg, 1993; Niebur and Erdös, 1991). They also show how muscle cells return proprioceptive signals back to the neuronal circuit simulator and how animals intentionally distribute body weight for locomotion patterns (Hu et al., 2009). Similarly, the kinematic simulator of fish (Ekeberg, 1993), which has a locomotion pattern in that the animal mostly undulates in a particular direction, was used with a neuronal network simulator to model the undulation of swimming behavior. This combination of kinematic simulator and neuronal network simulator was also used to model how the locomotion pattern changes due to a surrounding environment (Boyle et al., 2012) and how the central pattern generator arises from a few cells (Boyle et al., 2012; Izquierdo and Beer, 2018).

Even though there were studies on kinematic simulation of rod-shaped animals (Boyle et al., 2012; Ekeberg, 1993; Hu et al., 2009), to our best knowledge, there was no kinetic model that reproduces complex locomotion behavior of C. elegans, which includes all of the various modes of locomotion of C. elegans such as forward locomotion, backward locomotion, and turn from experimental observations. Instead, muscle cell activities from Ansatz (Hu et al., 2009), a hypothesis of the solution, or signals from a neuronal circuit simulator (Boyle et al., 2012; Ekeberg, 1993) were applied to the kinematic simulators. A simulator should have an operational structure that imitates physical quantities from an experiment to reproduce the motion of C. elegans in the experiment. However, until now, no kinetic simulation has such a structure. If there is a simulator that reproduces the motion of individual experiments, analysis of the kinetics of motion of specific experiments, which provides information on the individual force that exerts on each body part of the animal, will be enabled. Also, as the kinetic simulation reproduces the motion of C. elegans, the behavioral phenotype that emerged from the muscle activity of neuronal circuit simulation will be more credible.

We built a Newtonian-mechanics two-dimensional rigid body chain model of C. elegans to reproduce its locomotion. We incorporated its body angle, related to the contraction of the body wall muscle of C. elegans, into the primary operating principle of our kinetic model so that the model simulates measurable physical quantities of C. elegans from its experimental video. The model includes a chain of multiple rod rigid bodies, a damped torsional spring between the rigid bodies, and a control angle, which is the dynamic baseline angle from the value of the kymogram of a physical experiment. We formulated Newtonian equations of translational and rotational motion of the rigid body model and computed the numerical solution of the equations by numerical integration using the semi-implicit Euler method. As a result, we were able to demonstrate trajectories and kinetics of the general locomotion of C. elegans, such as crawling, swimming (Vidal-Gadea et al., 2011), omega-turn, and delta-turn (Broekmans et al., 2016).

Results

Newton’s equation of motion for locomotion of Caenorhabditis elegans: How does ElegansBot work?

We introduce the simple chain model of C. elegans' body. C. elegans has an elongated body along the head-to-tail axis. Thus, the worm’s body can be approximated as a midline extended along the anterial-posterial axis in the xy-coordinate plane (Figure 1A). Let $M$ (=2 µg, details in ‘Worm’s mass, actuator elasticity coefficient, and damping coefficient’ of Appendix) be the mass and $L$ (=1 mm) be the length of the worm. Midline was approximated as $n$ (=25) straight rods, whose ends are connected to the ends of neighboring rods (Figure 1A). The mass, length, and moment of inertia of each rod is $m = M / n$ , $2 r = L / n$ , and $I = m r^{2} / 3$ , respectively. When numbering the rods in order, with the rod at the end of the head being labeled as ‘1-rod’ and the rod at the end of the tail being labeled as ‘n-rod,’ let us designate the i-th rod as ‘i-rod.’ The point where i-rod and (i+1)-rod meets is ‘i-joint.’.

Figure 1

Download asset Open asset

Components of ElegansBot.

(A) Chain model for *C. elegans* body. (B) Rods in chain model. (C) i-actuator, which is a damped torsional spring. (D) Frictional force (black arrow) due to the translation motion of a rod. (E) Frictional force (black arrow) due to the rotational motion of a rod. (F) Joint force (black arrows) acting on i-rod and (i+1)-rod.

The motion of the worm corresponds to the motion of all the rods. To describe the motion of each rod (i-rod), we need to determine the displacement vector ( $d_{i}$ ), velocity vector ( $v_{i}$ ), the angle measured counterclockwise from the positive x-axis to the tangential direction of the rod ( $s_{i}$ ) (Figure 1B), and angular velocity ( $ω_{i}$ ) of i-rod at a given time $t$ . However, the minimum information required to describe the motion of all rods includes the displacement vector ( $d_{c}$ ) and velocity vector ( $v_{c}$ ) of the worm’s center of mass, $s_{i}$ and $ω_{i}$ for each rod (Details in ‘Minimum information required to describe the motion of each rod’ of Appendix).

Value of time-dependent variables such as $d_{c}$ , $v_{c}$ , $s_{i}$ , and $ω_{i}$ at a given time, $t$ will be expressed as $*^{(t)}$ . When initial values, $d_{c}^{(0)}$ , $v_{c}^{(0)}$ , $s_{i}^{(0)}$ , and $ω_{i}^{(0)}$ are given, the Newtonian equation of motion for acceleration, $a_{c}$ and angular acceleration, ${α_{i}}_{i \in 1, \dots, n}$ must be acquired and numerically integrated twice to find $d_{c}^{(t)}$ , $v_{c}^{(t)}$ , $s_{i}^{(t)}$ , and $ω_{i}^{(t)}$ at a given time, $t$ . To obtain the Newtonian equations of motion, we must find every force and torque acting on each rod. There are frictional force, muscle force, and joint force among types of forces acting on the rod, and there are frictional torque, muscle torque, and joint torque among types of torques whose descriptions are as follows.

The only external force acting on the worm is a frictional force from a ground surface such as an agar plate or water. The frictional force is an anisotropic Stokes frictional force, with a magnitude proportional to the speed and assumed different friction coefficients in perpendicular and parallel directions (Boyle et al., 2012), which guarantees that linearity in velocity is preserved in frictional force as well (Details in ‘Preservation of linearity in friction’ of Appendix). Because of this preservation of linearity, the frictional forces of translational motion (Figure 1D) and rotational motion (Figure 1E) can be calculated separately and added together to find total frictional force and torque. Previously known values of the friction coefficients in perpendicular and parallel directions are used (Boyle et al., 2012).

Let the perpendicular and parallel friction coefficients be $b_{⊥}$ and $b_{∥}$ for a straightened worm, respectively. Each rod experiences $1 / n$ of the frictional force the worm gets. Thus, the perpendicular and parallel friction coefficients of each rod are $b_{⊥} / n$ , $b_{∥} / n$ , respectively. The ratio of perpendicular friction coefficient to parallel friction coefficient ( $b_{⊥} / b_{∥}$ ) is 40 in agar plate and 1.5 in water (Berri et al., 2009; Boyle et al., 2012). This ratio is an important determining factor in whether the locomotion would be crawling or swimming (Boyle et al., 2012). The total frictional force that i-rod receives is $F_{b, i} = - \frac{b_{∥}}{n} (v_{i} \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - \frac{b_{⊥}}{n} (v_{i} \cdot {\hat{N}}_{i}) {\hat{N}}_{i}$ (‘·’: dot product of vectors, ${\hat{r}}_{i} \equiv {[\begin{matrix} \cos (s_{i}) & \sin (s_{i}) \end{matrix}]}^{T}$ : unit vector parallel to i-rod, ${\hat{N}}_{i} \equiv {[\begin{matrix} - \sin (s_{i}) & \cos (s_{i}) \end{matrix}]}^{T}$ : unit vector perpendicular to i-rod), and the total frictional torque that i-rod receives is $τ_{b, i} = - \frac{1}{3} \frac{b_{⊥}}{n} r^{2} ω_{i}$ (positive or negative values are for torque pointing away from or into the paper plane, respectively.) (The proof is in ‘Frictional torque by rotational motion’ of Appendix).

Mature hermaphrodite C. elegans has four muscle strands at the left dorsal, right dorsal, left ventral, and right ventral part of the body, and each muscle strand has 24, 24, 23, and 24 muscle cells, respectively (White et al., 1986). Muscle cells at similar positions on the anterior-posterior axis have an activity pattern in that muscles on one side (either dorsal or ventral) cooperate, and those on the opposite side have alternative activities. (Hall and Altun, 2008).

Therefore, we modeled a group of about four muscle cells, which are left dorsal, right dorsal, left ventral, and right ventral, at the same position on the anterior-posterior axis as one actuator (Figure 1C) so that there is a total of 24 ( $≃ (24 + 24 + 23 + 24) / 4$ ) actuators in the worm. On i-joint of the chain, there is an actuator labeled as i-actuator. As the number of actuators is 24, we set the number of rod( $n$ ) as 25, which is one more than the number of actuators. The actuator was modeled as a damped torsional spring due to the viscoelastic characteristics of muscle (Boyle et al., 2012; Hill, 1938). If the dorsal muscles of i-actuator contract more than the ventral muscles, i-actuator will bend to the dorsal direction and vice versa. To express this phenomenon by an equation, we defined the torque that i-actuator exerts on i-rod as $τ_{i} = τ_{κ, i} + τ_{c, i}$ that the elastic term is $τ_{κ, i} = κ (θ_{i} - θ_{ctrl, i})$ and the damping term is $τ_{c, i} = c (ω_{i + 1} - ω_{i})$ where $θ_{c t r l, i}$ is control angle, $θ_{i} = s_{i + 1} - s_{i}$ , and $κ$ and $c$ are the elasticity and damping coefficients of an actuator, respectively.

Control angle ( $θ_{c t r l, i}$ ) is a variable inside the elastic part of the muscle torque ( $τ_{κ, i}$ ), to which $τ_{κ, i}$ drives $θ_{i}$ close. Also, the control angle ( $θ_{c t r l, i}$ ), which can be expressed by a heatmap (Figures 2A, C, 3A and B), is an input value based on experimental data, a numerical model, or a neuronal network model. $τ_{c, i}$ represents the damping effect of muscle cells and somatic cells near i-actuator. The elasticity coefficient ( $κ$ ) and damping coefficient ( $c$ ) of an actuator were induced from previously known values (Boyle et al., 2012) (Details in ‘Worm’s mass, actuator elasticity coefficient, and damping coefficient’ of Appendix).

Figure 2 with 1 supplement see all

Download asset Open asset

Simulated locomotion from a sine kymogram.

(A) Crawling kymogram. Kymogram indicates the angle of i-joint which is located between i-rod and (i+1)-rod. Red and blue color mean i-joint bend in the dorsal and ventral directions, respectively. (B) Crawling trajectory. The yellow circle indicates the position of the worm’s head. The Orange and sky-blue lines show the worm’s head and tail trajectories, respectively. (C) Swimming kymogram. (D) Swimming trajectory.

Figure 3 with 2 supplements see all

Download asset Open asset

Simulated locomotion from a kymogram of a real worm locomotion video.

The length and direction of a black arrow indicate the magnitude and direction of the frictional force ( $- F_{b, i}^{(t)}$ ) that the corresponding body part, which is the starting point of the arrow, exerts on the surface. (A) Escaping behavior kymogram. Triangles over the heatmap indicate the corresponding time of snapshots shown in Figure (C). (B) Delta-turn kymogram. Triangles over the heatmap indicate the corresponding time of snapshots shown in Figure (D). (C) Escaping behavior trajectory. (D) Delta-turn trajectory. The arrow length scale is different from Figure (C) to clearly show the arrows' directions and head and tail tracks.

By assuming that i-rod receives torque ( $τ_{i}$ ) from i-actuator, and (i+1)-rod receives torque ( $- τ_{i}$ ), we can depict the bending that arises from the differential contraction of the dorsal and ventral muscles in i-actuator. The total muscle torque that i-rod receives from damped torsional springs on both ends is $τ_{c κ, i} = τ_{i} - τ_{i - 1}$ . The total muscle force ( $F_{c κ, i}$ ) that i-rod receives from both of its ends is as follows (Details in ‘Proof of muscle force’ of Appendix).

F_{c κ, i} = \sum_{j = - 1}^{0} \frac{{(- 1)}^{j - 1} τ_{i + j} \sin \frac{θ_{i + j}}{2}}{r \cos^{2} \frac{θ_{i + j}}{2}} [\begin{matrix} \cos (\frac{s_{i + j} + s_{i + j + 1}}{2}) \\ \sin (\frac{s_{i + j} + s_{i + j + 1}}{2}) \end{matrix}]

Two neighboring rods (i-rod and (i+1)-rod) are connected at i-joint. Therefore, when a force is applied to i-rod, (i+1)-rod also receives distributed force (Figure 1F) which we name as ‘joint force.’ The joint force that (i+1)-rod exerts on i-rod is symbolized as $F_{i} \equiv F_{(i + 1) i}$ . By Newton’s third law of motion about action and reaction, the joint force that i-rod exerts on (i+1)-rod is $F_{i (i + 1)} = - F_{(i + 1) i} = - F_{i}$ (Figure 1F). Joint force ( $F_{i}$ ) can be calculated from the previously introduced given values ( $s_{i}$ , $F_{c κ, i}$ , $F_{b, i}$ , $τ_{c κ, i}$ , $τ_{b, i}$ ) (Details in ‘Joint force calculation method’ of Appendix). When $F_{0} = F_{n} = 0$ , then the total joint force that i-rod receives is $F_{j o i n t, i} = F_{i} - F_{i - 1}$ and the total torque caused by joint force is $τ_{j o i n t, i} = [r_{i} \times (F_{i} + F_{i + 1})]$ where ‘×’ between two vectors means cross-product.

As all forces and torques are found, $d_{c}^{(t)}$ , $v_{c}^{(t)}$ , $s_{i}^{(t)}$ , $ω_{i}^{(t)}$ can be calculated by solving translational and rotational Newtonian equations of motion with numerical integration. The time-step ( $Δ t$ ) used in this work is $10^{- 5}$ s unless otherwise noted. Because the only external force exerts on the worm is the frictional force, the equation of translational motion is $a_{c} = \frac{\sum_{i} F_{b, i}}{M}$ . If friction coefficients ( $b_{⊥}$ , $b_{∥}$ ) are significantly greater than $\frac{M}{Δ t}$ , numerical integration using the explicit Euler method ( $v_{c}^{(t + Δ t)} = v_{c}^{(t)} + a_{c}^{(t)} Δ t = v_{c}^{(t)} + \frac{\sum_{i} F_{b, i}}{M} Δ t$ ) becomes unstable (Butcher, 2003). So, we tackled this instability of numerical integration by developing semi-implicit Euler Method ( $v_{c}^{(t + Δ t)} = v_{c}^{(t)} + a_{c}^{(t + Δ t)} Δ t ≃ v_{c}^{(t)} + \frac{1}{1 + \frac{b_{⊥} Δ t}{M}} \frac{\sum_{i} F_{b, i}^{(t)}}{M} Δ t$ ), which makes numerical integration stable when any frictional coefficients greater than or equal to 0 is given (Details in ‘Proof of numerical integration for the translational motion of a worm using semi-implicit Euler method’ of Appendix).

The equation of rotational motion of i-rod is $I α_{i} = τ_{t o t a l, i} = τ_{b, i} + τ_{c κ, i} + τ_{j o i n t, i}$ . When the friction-related value ( $b_{∥} r^{2}$ ), elasticity-related value ( $κ Δ t$ ), or damping coefficient( $c$ ) is significantly larger than $\frac{I}{Δ t}$ , numerical integration using explicit Euler method ( $ω_{i}^{(t + Δ t)} = ω_{i}^{(t)} + α_{i}^{(t)} Δ t = ω_{i}^{(t)} + \frac{τ_{t o t a l, i}^{(t)}}{I} Δ t$ ) becomes unstable (Butcher, 2003). To solve this instability, we constructed a semi-implicit Euler method for rotational motion and an error-corrected equation for angular momentum (Details in ‘Numerical integration of the rotational motion of i-rod using semi-implicit Euler method’ and ‘Correction formula for the rotational inertia of the entire worm’ of Appendix). By using these semi-implicit Euler methods, solutions for $d_{c}^{(t)}$ , $v_{c}^{(t)}$ , $s_{i}^{(t)}$ , $ω_{i}^{(t)}$ of a worm at a given time can be available for the ground surface of agar whose $b_{⊥}$ , $b_{∥}$ are significantly larger than $\frac{M}{Δ t}$ , water which has smaller friction coefficients than agar, or frictionless ground surface.

Can C. elegans in ElegansBot crawl or swim?

A kymogram is a heatmap that shows body angle, $θ_{i}^{(t)}$ ( $i \in {1, \dots, n - 1}$ ) at a given time, $t$ . By fitting a sine function to the kymogram of previous work (Vidal-Gadea et al., 2011), we obtained linear-wavenumber (after this referred to as wavenumber) and period of C. elegans crawling on the agar plate and swimming in water. The wavenumber ( $ν$ ) and the period ( $T$ ) are, respectively, 1.832 and 1.6 (s) on the agar plate and 0.667 and 0.4 (s) in water. For both crawling and swimming, amplitude ( $A$ ) was set to 0.6 (rad) arbitrarily to match the trajectory shown in the experimental video (Vidal-Gadea et al., 2011). Each kymogram of crawling (Figure 2A) and swimming (Figure 2C) was calculated by substituting amplitude ( $A$ ), wavenumber ( $ν$ ), and period ( $T$ ) into into $θ_{c t r l, i}^{(t)} = A \cos (2 π (ν (i - 1) / (n - 2) - t / T))$ .

Crawling trajectory, which performs sinusoidal locomotion in the positive x-axis direction, was obtained by inputting a crawling kymogram as $θ_{c t r l, i}^{(t)}$ input to ElegansBot (Figure 2B). Regarding crawling, the head track and the tail track have similar shapes. However, the tail track is more toward the negative x-axis direction than the head track. The difference between the head and tail tracks indicates that the worm pushes the ground surface by the distance between the head track and tail track to obtain thrust (Figure 2B). Indeed, we found that the body part placed diagonally with respect to the direction of the worm’s locomotion is pushing along the ground surface (Figure 2—figure supplement 1A). The thrust force of the worm cancels out most of the drag force, which enables the worm to move at nearly constant velocity. The average velocity of the worm is 0.208 (mm/s), which is consistent with the known values (Cohen et al., 2012; Jung et al., 2016; Omura et al., 2012; Shen et al., 2012).

In the previous work, the worm showed swimming behavior in a water droplet on an agar plate (Vidal-Gadea et al., 2011). As the friction coefficient of water is smaller than that of agar, even though the area that the worm swept was wider during swimming than crawling, the worm did not move forward much in comparison to the area it swept (Figure 2D). The worm gained significant momentum in the forward direction of locomotion when the body bent in the c-shape (Figure 2—figure supplement 1B). In contrast to crawling, during swimming, the worm did not receive constant thrust force over time. Thus, the speed of the worm exhibited significant oscillations over time (Figure 2—figure supplement 1B), and the average velocity was 0.223 (mm/s).

ElegansBot exhibits more complex behavior including the turn motion

Unlike previous C. elegans body kinematic simulation studies, our simulation can replicate the worm’s behavior using a kymogram (Figure 3A and B) derived from experimental videos. We utilized open-source software, Tierpsy Tracker (Javer et al., 2018) and WormPose (Hebert et al., 2021), to obtain the kymogram input ( $θ_{c t r l, i}^{(t)}$ ) for the ElegansBot. Through simulation, we aimed to reproduce the omega-turn and delta-turn behaviors observed in the experimental videos (Broekmans et al., 2016). When we used the vertical and horizontal friction coefficients $b_{⊥}$ and $b_{∥}$ on agar, as proposed in the previous work (Boyle et al., 2012), the trajectory was not accurately replicated. Given that the friction coefficients could vary depending on the concentration of the agar gel, we used $b_{⊥} / 100$ and $b_{∥} / 100$ for the vertical and horizontal friction coefficients, respectively, which resulted in a better trajectory replication (Details in ‘Proper selection of friction coefficients’ in Appendix).

The trajectory (Figure 3C and D, Figure 3—videos 1 and 2) obtained from ElegansBot accurately reproduces the experimental video (Broekmans et al., 2016). The changes in the direction of movement caused by turns are well replicated. Additionally, during the omega-turn or delta-turn, the body briefly performs a deep bend, and we newly discovered the mechanism that gains significant propulsion from the deep bend region to change direction using ElegansBot (Figure 3C and D). Moreover, the ElegansBot accurately reproduces not only the turns but also complex behaviors like the sequence of forward-backward-turn-forward, also known as escaping behavior.

Additionally, we calculated the mechanical power of the worm as a quantitative indicator to explain its locomotion during sequenced locomotive behavior, based on behavior classification (forward, backward locomotion, or turn, as defined in Methods). During escaping behavior, the worm produced an average power of 2094 fW in the initial forward locomotion, followed by an average of 16,437 fW (7.85 times that of the initial forward locomotion) in backward locomotion, and an average of 11,118 fW (5.31 times that of the initial forward locomotion) during turning (Figure 4A). After turning and resuming forward locomotion, it produced an average power of 5480 fW (2.62 times that of the initial forward locomotion). This indicates that the worm produced more power than that of initial forward locomotion to escape sudden threats. Let’s denote the average of a quantity for all given $i$ as ${⟨ * ⟩}_{i}$ . At the moment the worm formed a deep bend (t=11.2 s), the average magnitude of frictional force of the body part forming the deep bend ( ${⟨ | F_{b, i}^{(t)} | ⟩}_{i}$ where i=4 to 15) was 3536 pN, compared to the average magnitude of the remaining parts ( ${⟨ | F_{b, i}^{(t)} | ⟩}_{i}$ = 1737 pN where i=1 to 3 or i=16 to 25), which was 2.04 times greater (Figure 3C, Figure 4A). We analyzed delta-turn in the same manner. The worm produced an average power of 3514 fW in the initial forward locomotion, followed by an average of 11,176 fW (3.18 times that of the initial forward locomotion) in subsequent backward locomotion. In the relatively short duration of forward locomotion following the backward locomotion, the worm produced an average power of 17,544 fW (4.99 times that of the initial forward locomotion), and an average of 13,046 fW (3.71 times that of the initial forward locomotion) during turns (Figure 4B). After the turn, when resuming forward locomotion, the worm produced an average power of 6429 fW (1.83 times that of the initial forward locomotion). At the moment the worm formed a deep bend (t=6.1 s), the average magnitude of frictional force of the body part ( ${⟨ | F_{b, i}^{(t)} | ⟩}_{i}$ where i=16 to 25) was 10,497 pN, compared to the average magnitude of the remaining parts ( ${⟨ | F_{b, i}^{(t)} | ⟩}_{i}$ = 2,677 pN where i=1 to 15), which was 3.92 times greater (Figure 3D, Figure 4B). In both escaping behavior and delta-turn, the worm consistently produced more power in the subsequent backward locomotion and turn than in the initial forward locomotion.

Figure 4 with 1 supplement see all

Download asset Open asset

Frictional force on each rod.

(A) Escaping behavior. The top panel represents the frictional force $F_{b, i}^{(t)}$ experienced by i-rod. As indicated on the color wheel to the right, the hue of this heatmap represents the direction of the force, and the saturation represents the magnitude of the force. The second panel from the top shows the magnitude of the frictional force $| F_{b, i}^{(t)} |$ . The third panel from the top represents the average ${⟨ | F_{b, i}^{(t)} | ⟩}_{i}$ (black solid line) of each column in the middle panel and the power (red dotted line), which is the amount of energy the worm consumes per unit time. The bottom panel represents the classification of the worm’s behavior (blue: forward locomotion, red: backward locomotion, green: turn) (The definitions of behavioral categories are in Methods). The triangles over each panel indicate the corresponding time of the snapshots depicted in Figure 3C. (B) Delta-turn. The triangles over each panel indicate the corresponding time of the snapshots depicted in Figure 3D.

ElegansBot presents body shape ensembles of C. elegans from a shape in water en route to agar

While there have been studies on how locomotion patterns change in agar and water by merging neural and kinematic simulations (Boyle et al., 2012), there have been none that solely used kinetic simulation to analyze how speed manifests depending on the frequency and period of locomotion. We demonstrate this aspect. We studied the locomotion speed of the worm under different friction coefficients, which represent the influence of water, agar, and intermediate frictional environment, using ElegansBot. The vertical and horizontal friction coefficients in water are $b_{w a t e r, ⊥} = 5.2 \times 10^{3}$ (μg/sec) and $b_{w a t e r, ∥} = b_{w a t e r, ⊥} / 1.5$ , respectively, while in agar, these values are $b_{a g a r, ⊥} = 1.28 \times 10^{8}$ (μg/sec) and $b_{a g a r, ∥} = b_{a g a r, ⊥} / 40$ (Boyle et al., 2012). For environmental index $σ \in [0,1]$ , we have defined the vertical and horizontal friction coefficients in the environment between water ( $σ = 0$ ) and agar ( $σ = 1$ ) as $b_{σ, ⊥} = b_{w a t e r, ⊥}^{1 - σ} b_{a g a r, ⊥}^{σ}$ and $b_{σ, ∥} = b_{w a t e r, ∥}^{1 - σ} b_{a g a r, ∥}^{σ}$ , respectively.

Under an environmental index $σ$ , for various pairs of frequency-period ( $ν$ , $T$ ) when the control angle is $θ_{c t r l, i}^{(t)} = A \cos (2 π (ν \frac{i - 1}{23} - \frac{t}{T}))$ (with $A$ = 0.6 (rad)), we have found the ( $ν$ , $T$ ) that maximizes the worm’s average velocity(optimal ( $ν$ , $T$ )) (Figure 5A). The optimal ( $ν$ , $T$ ) exhibits a nearly linear distribution (Figure 5B). We noticed a transition from swimming body shape to crawling body shape as $σ$ varies (Figure 5, Figure 5—figure supplement 1). The optimal ( $ν$ , $T$ ) for $σ$ =0(water) is (0.65, 0.4 s), matching the actual ( $ν$ , $T$ ) value of swimming behavior (Vidal-Gadea et al., 2011). The optimal ( $ν$ , $T$ ) for $σ$ =1(agar) is (1.9, 0.8 s), and the optimal $ν$ (1.9) matches the actual $ν$ value (1.832) for crawling behavior (Vidal-Gadea et al., 2011), with the optimal T (0.8 s) being half the actual T value (1.6 s).

Figure 5 with 3 supplements see all

Download asset Open asset

Body shape transition from the shape in water to the shape in agar.

(A) Average velocity of the worm as a function of wavenumber ( $ν$ ) and period ( $T$ ) for a given friction coefficient. The star symbol indicates the pair of ( $ν$ , $T$ ) that maximizes the worm’s average velocity. (B) For each environmental index $σ$ , the pair of ( $ν$ , $T$ ) that maximizes the worm’s average velocity. The worm figures inside the small rectangles pointed to by the arrows represent the body shape corresponding to the respective ( $ν$ , $T$ ) pair.

We wanted to understand the impact of the environmental index $σ$ not only on forward locomotion but also on sequenced locomotive behavior. First, we analyzed the effect of the environmental index $σ$ on escaping behavior as follows. Let’s denote the set of a quantity for all pairs of index $i$ and time $t$ as ${*}_{i, t}$ . When the escaping behavior kymogram input ${θ_{c t r l, i}^{(t)}}_{i, t}$ was same as Figure 3A, we explored the effect of vertical and horizontal friction coefficients on the worm’s motion. Where $σ$ ranged from 1.0 to 0, the trajectory varied with $σ$ (Figure 5—figure supplement 2A), and $E_{θ} = {⟨ | θ_{i}^{(t)} - θ_{c t r l, i}^{(t)} | ⟩}_{i, t}$ decreased as $σ$ decreased (Figure 5—figure supplement 2B). From $σ$ = 1.0 to $σ$ = 0.1, the total absolute angular change ( $S = \sum_{t = 0}^{T - Δ t} | {⟨ s_{i}^{(t + Δ t)} ⟩}_{i} - {⟨ s_{i}^{(t)} ⟩}_{i} |$ where $T$ is the total time of the experimental video.) increased as $σ$ decreased. However, from $σ$ = 0.7 to $σ$ = 0, $S$ remained constant within the error of 0.33 rad, and the total traveled distance( $\sum_{t} | v_{c}^{(t)} Δ t |$ ) decreased as $σ$ decreased. The maximum total traveled distance was at $σ$ = 0.8. Using the same analysis method with the kymogram input ${θ_{c t r l, i}^{(t)}}_{i, t}$ same as Figure 3B, we analyzed the impact of the environmental index $σ$ on delta-turn. Where $σ$ ranged from 1.0 to 0, the trajectory varied with $σ$ (Figure 5—figure supplement 3A), and $E_{θ}$ also decreased as $σ$ decreased (Figure 5—figure supplement 3B). From $σ$ = 1 to $σ$ = 0.6, $S$ increased as $σ$ decreased. From $σ$ = 0.6 to $σ$ = 0, $S$ decreased as $σ$ decreased. The maximum total traveled distance was at $σ$ = 0.9. From $σ$ = 0.9 to $σ$ = 0, the total traveled distance decreased as $σ$ decreased.

Discussion

ElegansBot is an advanced kinetic simulator that reproduces C. elegans' various locomotion

The known crawling speed range of C. elegans (Cohen et al., 2012; Jung et al., 2016; Omura et al., 2012; Shen et al., 2012) matches the speed in our simulation. The force dispersion pattern of the forward movement of a snake (Hu et al., 2009) is similar to the force dispersion pattern of crawling in our model, where the body part placed diagonally in the direction of movement generates thrust. The head and tail tracks of our simulation resemble the trace left on the agar plate by C. elegans during locomotion (Yeon et al., 2018), providing evidence of the mechanism where C. elegans moves forward by pushing along the ground surface. Given that friction and elasticity coefficients can vary between experiments, the appropriate selection of these values allows the trajectories of omega-turns and delta-turns in our simulations to match the experimental videos (Broekmans et al., 2016). Previous work (Berri et al., 2009; Boyle et al., 2012) eliminated inertia from the equations of motion, but our simulation includes it, allowing calculation even in cases where inertia is significant due to low friction coefficients. Using the crawling and swimming wavenumbers and periods from the experiments (Vidal-Gadea et al., 2011), we computed sine functions to create trajectories for crawling and swimming. We also analyzed how friction forces act on the worm during crawling and swimming, studying how the worm gains propulsion. We demonstrated that we could reproduce various locomotion observed in experimental videos, such as forward-backward-(omega turn)-forward constituting escaping behavior and delta-turn navigation, by providing the kymogram obtained from representative physical values from the experimental videos, as well as the kymogram obtained from a program (Hebert et al., 2021; Javer et al., 2018) extracting the body angles from actual experimental videos into ElegansBot. Our established Newtonian equations of motion are accurate and robust, suggesting that not only does our simulation replicate the experimental videos, but it also provides credible estimates for detailed forces.

ElegansBot will serve as a strong bridge for enhancing the knowledge in ‘from-synapse-to-behavior’ research

Our method could be used for kinetic analysis of behaviors not covered in this paper. It could also be used when analyzing behavior changes caused by mutation or ablation experiments. Given that our simulation allows for kinetic analysis, it could be used to calculate the energy expended by the worm during locomotion, serving as an activity index. Our simulation only requires the body angles as input data, so even if the video angle shifts and trajectory information is lost, the trajectory can be recovered from the kymogram. Our simulation could also be used when studying neural circuit models of C. elegans. It could be used to check how signals from neural network models manifest as behaviors which is a needed function from previous work (Sakamoto et al., 2021), and it could be used when studying compound models of neural circuits and bodies. For example, when creating models that receive proprioception input based on body shape (Boyle et al., 2012; Ekeberg, 1993; Izquierdo and Beer, 2018; Niebur and Erdös, 1991), our method could be used. Finally, our method could be used in general for the broad utility to analyze the motion of rod-shaped animals like snakes or eels and to simulate the motion of rod-shaped robots.

Methods

Frequency and wavelength of C. elegans locomotion

Sine function fitting was applied to the crawling and swimming kymograms (Vidal-Gadea et al., 2011) to determine the frequency and wavelength of C. elegans locomotion on agar and water.

C. elegans locomotion videos

Videos of C. elegans' escaping behavior and foraging behavior were obtained from previous work (Broekmans et al., 2016). A single representative video out of a total of one hundred escaping behavior videos was used as data in this paper. Additionally, only the delta-turning portion of the foraging behavior videos was cut out and used as data in this paper.

Obtaining kymograms from video

The following method was used to extract the kymogram from the video of C. elegans: The body angles and midline were extracted using Tierpsy Tracker (Javer et al., 2018) from the original video where the worm is locomoting. Tierpsy Tracker failed to extract the midline of the worm when the body parts meet or the worm is coiled. The midline information of the frames successfully predicted by Tierpsy Tracker and the original video information were used as ground truth training data for a program called WormPose (Hebert et al., 2021). WormPose trained an artificial neural network to extract the body angles and midline of a coiled worm using a generative method based on the input data. The body angles were extracted from the original video using the trained WormPose program. For the frames where Tierpsy Tracker failed to extract the body angles, it was replaced with the body angles extracted by WormPose. Nonetheless, there were frames where the body angle extraction failed. If the period of failed body angle prediction was continuously less than three frames (about 0.01 s), the body angles for that period was predicted using linear interpolation.

Program code and programming libraries

Equations for the chain model, friction model, muscle model, and numerical integration that constitute the ElegansBot were designed from the body angle information and kymogram that change every moment of time. Python (van Aken et al., 1995) version 3.8 was used to implement the equations constituting ElegansBot as a program. NumPy (Harris et al., 2020) version 1.19 was used for numerical calculations, and Numba (Lam et al., 2015) version 0.54 was used for CPU calculation acceleration. SciPy (Virtanen et al., 2020) version 1.5 was used for curve fitting and Savitzky-Golay filter (Savitzky and Golay, 1964) to classify the worm’s behavioral categories. The Matplotlib (Hunter, 2007) library was used to represent C. elegans' body pose and trajectory in figures and videos. The program code used in the research can be obtained from the open database GitHub (Taegon Chung, 2023, ElegansBot, 1.0.1, https://github.com/taegonchung/elegansbot, copy archived at Chung, 2024) or Python software repository PyPI (https://pypi.org/project/ElegansBot/), and the web live demo can be found at GitHub Page (https://taegonchung.github.io/elegansbot/). This code calculated the ElegansBot simulation of 10 s of simulation time in about 10 s of run-time on an Intel E3-1230v5 CPU.

Physical constants of the ground surface

The friction coefficient values for the ground surface where C. elegans crawled and swam and the elastic and damping coefficients of C. elegans muscles were obtained from previous work (Boyle et al., 2012). The muscle elasticity and damping coefficients were converted into coefficients for the damped torsional spring to be used in our model (Details in ‘Worm’s mass, actuator elasticity coefficient, and damping coefficient’ of Appendix).

Defining behavioral categories

We determined the classification of the worm’s behavior over time as follows. Let $ξ_{i} \equiv (i - 1) / (n - 2)$ (i.e., $0 \leq ξ_{i} \leq 1$ ). We defined a sine function fitting for the body angle $θ_{i}^{(t)}$ as ${\hat{θ}}_{i}^{(t)} = A^{(t)} \sin ((2 π / λ^{(t)}) (ξ_{i} - ξ_{0}^{(t)})) + θ_{0}^{(t)}$ (where $A^{(t)} \geq 0$ and $0 \leq ξ_{0}^{(t)} < λ^{(t)}$ ). Let us denote the set of a quantity for all $i$ as ${*}_{i}$ . For a given time $t$ , by curve fitting the function ${\hat{θ}}_{i}^{(t)}$ to the set of body angles ${θ_{i}^{(t)}}_{i}$ , we can obtain the parameters ( $A^{(t)}$ , $λ^{(t)}$ , $ξ_{0}^{(t)}$ , $θ_{0}^{(t)}$ ) (Figure 4—figure supplement 1). For curve fitting, we used 'curve_fit' from the scipy library (Virtanen et al., 2020). ‘curve_fit’ requires the function to be fitted, the data, and the initial guess values of the function parameters. Therefore, we determined the initial guess values ( ${\hat{A}}^{(t)}$ , ${\hat{λ}}^{(t)}$ , ${\hat{ξ}}_{0}^{(t)}$ , ${\hat{θ}}_{0}^{(t)}$ ) as follows. For $t = 0$ , we calculated ( ${\hat{A}}^{(0)}$ , ${\hat{λ}}^{(0)}$ , ${\hat{ξ}}_{0}^{(0)}$ , ${\hat{θ}}_{0}^{(0)}$ ) using the following equations:

\begin{array}{ll} {\hat{A}}^{(0)} = & \frac{1}{2} (max_{i} (θ_{i}^{(0)}) - min_{i} (θ_{i}^{(0)})) \\ {\hat{λ}}^{(0)} = & \frac{2}{n - 2} (\underset{i}{a r g m a x} (θ_{i}^{(0)}) - \underset{i}{a r g m i n} (θ_{i}^{(0)})) \\ {\hat{ξ}}_{0}^{(0)} = & \frac{1}{n - 2} (\underset{i}{a r g m a x} (θ_{i}^{(0)})) - \frac{{\hat{λ}}^{(0)}}{4} \\ {\hat{θ}}_{0}^{(0)} = & \frac{1}{n - 1} (\sum_{i} θ_{i}^{(0)}) \end{array}

Using these initial guess values, we curve-fitted ${\hat{θ}}_{i}^{(0)}$ to ${θ_{i}^{(0)}}_{i}$ to obtain ( $A^{(0)}$ , $λ^{(0)}$ , $ξ_{0}^{(0)}$ , $θ_{0}^{(0)}$ ). For $t \geq Δ t$ , we obtained the initial guess values for curve fitting as ${\hat{A}}^{(t)} = A^{(t - Δ t)}$ , ${\hat{λ}}^{(t)} = λ^{(t - Δ t)}$ , ${\hat{ξ}}_{0}^{(t)} = ξ_{0}^{(t - Δ t)}$ , ${\hat{θ}}_{0}^{(t)} = θ_{0}^{(t - Δ t)}$ . Then, using these initial guess values, we curve-fitted ${\hat{θ}}_{i}^{(t)}$ to ${θ_{i}^{(t)}}_{i}$ to obtain ( $A^{(t)}$ , $λ^{(t)}$ , $ξ_{0}^{(t)}$ , $θ_{0}^{(t)}$ ). Since the phase $ξ_{0}^{(t)}$ is not continuous for all time $t$ , we defined a continuous value ${\tilde{ξ}}_{0}^{(t)}$ for all time $t$ as follows (Figure 4—figure supplement 1):

{\tilde{ξ}}_{0}^{(t)} = {\begin{cases} ξ_{0}^{(t)} & i f t = 0 \\ ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} + {\tilde{ξ}}_{0}^{(t - Δ t)} & if t > 0 a n d - \frac{λ^{(t)}}{2} \leq ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} \leq \frac{λ^{(t)}}{2} \\ ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} + {\tilde{ξ}}_{0}^{(t - Δ t)} + λ^{(t)} & if t > 0 a n d ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} < - \frac{λ^{(t)}}{2} \\ ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} + {\tilde{ξ}}_{0}^{(t - Δ t)} - λ^{(t)} & if t > 0 a n d \frac{λ^{(t)}}{2} < ξ_{0}^{(t)} - ξ_{0}^{(t - Δ t)} \end{cases}

To obtain the derivatives of the noise-reduced smoothed values ${\bar{ξ}}_{0}^{(t)}$ and ${\bar{θ}}_{0}^{(t)}$ for the raw data ${\tilde{ξ}}_{0}^{(t)}$ and $θ_{0}^{(t)}$ , respectively, we applied a Savitzky-Golay filter (Savitzky and Golay, 1964; Virtanen et al., 2020). This filter, set with a smoothing time window of 0.5 s, a polynomial order of 2, and a derivative order of 1, yielded $\frac{d {\bar{ξ}}_{0}^{(t)}}{d t}$ and $\frac{d {\bar{θ}}_{0}^{(t)}}{d t}$ . We then calculated the temporal integrals ${ξ^{'}}_{0}^{(ζ)} \equiv \sum_{t = 0}^{ζ} \frac{d {\bar{ξ}}_{0}^{(t)}}{d t} Δ t$ and ${θ^{'}}_{0}^{(ζ)} \equiv \sum_{t = 0}^{ζ} \frac{d {\bar{θ}}_{0}^{(t)}}{d t} Δ t$ . Let us denote the average of a quantity for all time t as ${⟨ * ⟩}_{t}$ . We calculated ${\bar{ξ}}_{0}^{(t)} = {ξ^{'}}_{0}^{(t)} - {⟨ {ξ^{'}}_{0}^{(t)} ⟩}_{t} + {⟨ {\tilde{ξ}}_{0}^{(t)} ⟩}_{t}$ and ${\bar{θ}}_{0}^{(t)} = {θ^{'}}_{0}^{(t)} - {⟨ {θ^{'}}_{0}^{(t)} ⟩}_{t} + {⟨ θ_{0}^{(t)} ⟩}_{t}$ (Figure 4—figure supplement 1). Finally, we defined the worm’s behavior classification as turn when ${\bar{θ}}_{0}^{(t)} < - 0.07$ , as forward locomotion when ${\bar{θ}}_{0}^{(t)} \geq - 0.07$ and $\frac{d {\bar{ξ}}_{0}^{(t)}}{d t} > 0$ , and as backward locomotion when ${\bar{θ}}_{0}^{(t)} \geq - 0.07$ and $\frac{d {\bar{ξ}}_{0}^{(t)}}{d t} \leq 0$ (Figure 4—figure supplement 1).

Appendix 1

Worm’s mass, actuator elasticity coefficient, and damping coefficient

1. Mass

The maximum radius of the worm is 40µm (Boyle et al., 2012). Assuming that the border surrounding the cross-section parallel to the anterior-posterior axis is a sine function, the average radius is $γ = 40 \times 2 / π ≃ 25 μ m$ . Assuming that the worm is a cylindrical body with a bottom surface radius of 25µm, the volume of the worm is $1 m m {(0.025 m m)}^{2} π ≃ 0.002 {m m}^{3}$ , and the density of the worm is $≃ 1000 μ g / m m^{3}$ (Reina et al., 2013), so the weight of the worm is $M = 1000 μ g / {m m}^{3} \times 0.002 {m m}^{3} = 2 μ g$ .

2. Torque elasticity coefficient

In previous research, the muscle elasticity and damping coefficient were designed as functions of the input signal (Boyle et al., 2012). The maximum value of this muscle elasticity coefficient is $k_{m a x} = 2.8 \cdot 10^{8} [μ g / {s e c}^{2}]$ , and the maximum value of the muscle damping coefficient is $(k_{m a x} / 5.6) \cdot (1 s e c)$ . When $θ_{i - 1} = θ_{i} = θ_{i + 1} = 0$ and the length of the moment arm where the muscle exerts force is equal to the average radius $γ$ of the worm, the change in the elastic torque due to the change in $θ_{i}$ is the torque elasticity coefficient $κ$ , so $κ = \frac{d τ_{κ, i}}{d θ_{i}} = \frac{d}{d θ_{i}} (γ (k (2 γ t a n (θ_{i} / 2)))) ≃ k γ^{2} = 1.75 \cdot 10^{5} μ g \cdot m m^{2} / (s e c^{2} \cdot rad)$ . In ElegansBot, it was assumed that this $κ$ value is constant regardless of $θ_{i - 1}$ , $θ_{i}$ , $θ_{i + 1}$ . In the same way, the torque muscle damping coefficient is $c = (1 / 5.6) \cdot 1.75 \cdot 10^{5} μ g \cdot {mm}^{2} / (s e c \cdot rad)$ .

Minimum information required to describe the motion of each rod

1. Information required to describe the movement of the rod

To describe the motion of all rods, it is necessary to know the position ( $d_{i}$ ), velocity ( $v_{i}$ ), angle ( $s_{i}$ ) measured counterclockwise from the positive x-axis to the direction of vector $r_{i}$ , and angular velocity( $ω_{i}$ ) of every i-rod. However, knowing only the position and velocity of the worm ( $d_{c}$ , $v_{c}$ ) and the angle and angular velocity of every i-rod ( $s_{i}$ , $ω_{i}$ ) is sufficient to describe the motion of all rods.

2. Boundary conditions for the position of the rod

The center of mass of the worm is $d_{c} = {[\begin{matrix} x_{c} & y_{c} \end{matrix}]}^{T} = (\sum_{i = 1}^{n} m d_{i}) / M = (\sum_{i = 1}^{n} d_{i}) / n$ , and the vector parallel to i-rod is $r_{i} \equiv r {[\begin{matrix} \cos (s_{i}) & \sin (s_{i}) \end{matrix}]}^{T} = r {\hat{r}}_{i}$ . If the position vector where i-rod and (i+1)-rod meet is $d_{{tip}_{i}}$ , and the free end of 1-rod and n-rod are $d_{{tip}_{0}}$ and $d_{{tip}_{n}}$ , then $d_{{tip}_{i}} = d_{i} + r_{i} = d_{i + 1} - r_{i + 1}$ (Appendix 1—figure 1A).

3. Method of calculating the relative and absolute position of the rod

If the relative coordinate with $d_{{tip}_{0}}$ as the origin is designated as $d `$ , then the following equations satisfy.

\begin{array}{ll} d_{{tip}_{0}}^{'} = [\begin{matrix} 0 \\ 0 \end{matrix}] \\ d_{{tip}_{i}}^{'} = \sum_{j = 1}^{i} 2 r_{j} \\ d_{i}^{'} = \frac{d_{{tip}_{i}}^{'} + d_{{tip}_{i - 1}}^{'}}{2} \\ d_{c}^{'} = \frac{1}{n} \sum_{i = 1}^{n} d_{i}^{'} \\ d_{i} = d_{i}^{'} - d_{c}^{'} + d_{c} \end{array}

Thus, $d_{i}$ can be calculated from $d_{c}$ and $s_{i}$ (Appendix 1—figure 1B).

4. Method of calculating the relative and absolute velocity of the rod

Based on the relationship between the worm’s momentum and the rods' momentum, $v_{c} = (\sum_{i = 1}^{n} m v_{i}) / M = (\sum_{i = 1}^{n} v_{i}) / n$ . In the same way as the location information compression, the relative velocity vector $v `$ with $v_{{tip}_{0}}$ as the origin has the following relationships ( $\hat{z} \equiv {\hat{r}}_{i} \times {\hat{N}}_{i}$ ).

\begin{array}{ll} v_{{tip}_{0}}^{'} & = [\begin{matrix} 0 \\ 0 \end{matrix}] \\ v_{{tip}_{i}}^{'} & = \sum_{j = 1}^{i} ω_{j} \hat{z} \times (2 r_{j}) \\ v_{i}^{'} & = \frac{v_{{tip}_{i}}^{'} + v_{{tip}_{i - 1}}^{'}}{2} \\ v_{c}^{'} & = \frac{1}{n} \sum_{i = 1}^{n} v_{i}^{'} \\ v_{i} & = v_{i}^{'} - v_{c}^{'} + v_{c} \end{array}

Thus, $v_{i}$ can be calculated from $v_{c}$ , $s_{i}$ , and $ω_{i}$ (Appendix 1—figure 1C).

Preservation of linearity in friction

The velocity $v$ of an arbitrary point particle which is included by i-rod can be decomposed into two velocity components $v = v_{α} + v_{β}$ . In this case, if the object is subject to an anisotropic Stokes friction, there is linearity between the friction $F_{b, α}$ , $F_{b, β}$ obtained from each velocity component $v_{α}$ , $v_{β}$ and the friction $F_{b}$ obtained from velocity $v$ .

\begin{array}{lll} F_{b, α} & = - b_{∥} n^{- 1} (v_{α} \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} (v_{α} \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ F_{b, β} & = - b_{∥} n^{- 1} (v_{β} \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} (v_{β} \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ F_{b} & = - b_{∥} n^{- 1} (v \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} (v \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ = - b_{∥} n^{- 1} ((v_{α} + v_{β}) \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} ((v_{α} + v_{β}) \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ = - b_{∥} n^{- 1} (v_{α} \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} (v_{α} \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ - b_{∥} n^{- 1} (v_{β} \cdot {\hat{r}}_{i}) {\hat{r}}_{i} - b_{⊥} n^{- 1} (v_{β} \cdot {\hat{N}}_{i}) {\hat{N}}_{i} \\ = F_{b, α} + F_{b, β} \end{array}

Frictional torque by rotational motion

Let us denote variable $ρ$ as the distance from the center of i-rod measured along the direction of vector ${\hat{r}}_{i}$ . It is to be noted that $ρ$ is within the range $[- r, r]$ . For an infinitesimal $d ρ$ where $0 < d ρ ≪ 1$ , the moment arm vector for the infinitesimal interval $[ρ - d ρ / 2, ρ + d ρ / 2]$ (hereafter referred to as the infinitesimal interval $ρ$ ) from the center of i-rod is $ρ {\hat{r}}_{i}$ . The coefficient of friction for the infinitesimal interval $ρ$ is:

\frac{b_{⊥}}{n} \frac{d ρ}{2 r}

The velocity component due to the rotational motion of the infinitesimal interval $ρ$ is $ρ ω_{i} {\hat{N}}_{i}$ . Therefore, the frictional force received by the infinitesimal interval $ρ$ due to rotational motion is:

- \frac{b_{⊥}}{n} \frac{d ρ}{2 r} ρ ω_{i} {\hat{N}}_{i} = - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} ρ d ρ {\hat{N}}_{i}

The torque received by the infinitesimal interval $ρ$ due to rotational motion is:

\begin{array}{ll} d τ & = (ρ {\hat{r}}_{i}) \times (- \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} ρ d ρ {\hat{N}}_{i}) \\ = - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} ρ^{2} d ρ {\hat{r}}_{i} \times {\hat{N}}_{i} \\ = - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} ρ^{2} d ρ \hat{z} \end{array}

The total frictional torque received by i-rod due to rotational motion is

\begin{array}{ll} τ_{b, i} & = \int_{- r}^{r} d τ \\ = \int_{- r}^{r} - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} ρ^{2} d ρ \hat{z} \\ = - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} {[\frac{1}{3} ρ^{3}]}_{- r}^{r} \hat{z} \\ = - \frac{1}{2} \frac{b_{⊥}}{n r} ω_{i} \frac{2}{3} r^{3} \hat{z} \\ = - \frac{1}{3} \frac{b_{⊥}}{n} r^{2} ω_{i} \hat{z} \end{array}

Proof of muscle force

The muscle force, which makes the torque $τ_{i}$ received from i-actuator to i-rod, was designed as follows. The damped torsion spring is connected at the center of each rod and gives a force in a direction perpendicular to the rod. The support is connected to i-actuator’s midpoint and i-joint (Appendix 1—figure 2A). Let us say that the mass of the support and i-actuator are both 0. The support gives forces ( $F_{{sup}_{1}, i}$ , $F_{{sup}_{2}, i}$ ) of the same size to i-rod and (i+1)-rod in a direction parallel to the support (Appendix 1—figure 2B). Let us assume that the resultant force received by i-actuator is $0$ (Appendix 1—figure 2C). Thus, the forces from the actuator generate no torque at the connection point in the middle of the rod but at the connection point at the end (Appendix 1—figure 2D).

\begin{array}{ll} F_{{sup}_{1}, i} & = F_{p_{1}, i} \cos \frac{θ_{i}}{2} \\ τ_{i} & = (r_{i} \times F_{{sup}_{1}, i}) \cdot \hat{z} \\ = r F_{{sup}_{1}, i} \cos \frac{θ_{i}}{2} \\ F_{{sup}_{1}, i} & = \frac{τ_{i}}{r \cos \frac{θ_{i}}{2}} \\ F_{p_{1}, i} & = \frac{F_{{sup}_{1}, i}}{\cos \frac{θ_{i}}{2}} \\ = \frac{τ_{i}}{r \cos^{2} \frac{θ_{i}}{2}} \end{array}

The force that i-rod receives from i-actuator is $F_{p_{1}, i} + F_{{sup}_{1}, i}$ and the magnitude of the force is $|F_{p_{1}, i} + F_{{sup}_{1}, i}| = F_{p_{1}, i} s i n (θ_{i} / 2) = τ_{i} s i n (θ_{i} / 2) / (r {c o s}^{2} (θ_{i} / 2))$ and the direction of the force is $- {[\begin{matrix} \cos ((s_{i} + s_{i + 1}) / 2) & \sin ((s_{i} + s_{i + 1}) / 2) \end{matrix}]}^{T}$ . These equations come down to the following equations.

F_{p_{1}, i} + F_{{sup}_{1}, i} = - \frac{τ_{i} \sin \frac{θ_{i}}{2}}{r \cos^{2} \frac{θ_{i}}{2}} [\begin{matrix} \cos (\frac{s_{i} + s_{i + 1}}{2}) \\ \sin (\frac{s_{i} + s_{i + 1}}{2}) \end{matrix}]

F_{p_{1}, i} + F_{{sup}_{1}, i} + F_{p_{2}, i} + F_{{sup}_{2}, i} = 0

F_{p_{2}, i} + F_{{sup}_{2}, i} = \frac{τ_{i} \sin \frac{θ_{i}}{2}}{r \cos^{2} \frac{θ_{i}}{2}} [\begin{matrix} \cos (\frac{s_{i} + s_{i + 1}}{2}) \\ \sin (\frac{s_{i} + s_{i + 1}}{2}) \end{matrix}]

Therefore, the total force that i-rod receives from i-actuator and (i-1)-actuator is

F_{c κ, i} = \sum_{j = - 1}^{0} \frac{{(- 1)}^{j - 1} τ_{i + j} \sin \frac{θ_{i + j}}{2}}{r \cos^{2} \frac{θ_{i + j}}{2}} [\begin{matrix} \cos (\frac{s_{i + j} + s_{i + j + 1}}{2}) \\ \sin (\frac{s_{i + j} + s_{i + j + 1}}{2}) \end{matrix}]

Joint force calculation method

Let us calculate the joint force $F_{i}$ from the given values ( $s_{i}$ , $F_{c κ, i}$ , $F_{b, i}$ , $τ_{c κ, i}$ , $τ_{b, i}$ ). In this part, the superscript $*^{T}$ means the transpose of a vector or a matrix. i-rod and (i+1)-rod always meet at i-joint can be expressed as a following vector equation.

d_{i} + r_{i} = d_{i + 1} - r_{i + 1}

Differentiating this equation twice for time, we can see that the accelerations of the ends of i-rod and (i+1)-rod at i-joint are the same.

\begin{array}{ll} \frac{d^{2}}{d t^{2}} (d_{i} + r_{i}) & = \frac{d^{2}}{d t^{2}} (d_{i} + r {\hat{r}}_{i}) \\ = \frac{d}{d t} (v_{i} + r ω_{i} {\hat{N}}_{i}) \\ = a_{i} + r \frac{d ω_{i}}{d t} {\hat{N}}_{i} + r ω_{i} \frac{d {\hat{N}}_{i}}{d t} \\ = a_{i} + r α_{i} {\hat{N}}_{i} - r ω_{i}^{2} {\hat{r}}_{i} \\ = a_{i} + α_{i} \times r_{i} - ω_{i}^{2} r_{i} \end{array}

\begin{array}{ll} \frac{d^{2}}{d t^{2}} (d_{i + 1} - r_{i + 1}) & = \frac{d^{2}}{d t^{2}} (d_{i + 1} - r {\hat{r}}_{i + 1}) \\ = \frac{d}{d t} (v_{i + 1} - r ω_{i + 1} {\hat{N}}_{i + 1}) \\ = a_{i + 1} - r \frac{d ω_{i + 1}}{d t} {\hat{N}}_{i + 1} - r ω_{i + 1} \frac{d {\hat{N}}_{i + 1}}{d t} \\ = a_{i + 1} - r α_{i + 1} {\hat{N}}_{i + 1} + r ω_{i + 1}^{2} {\hat{r}}_{i + 1} \\ = a_{i + 1} - α_{i + 1} \times r_{i + 1} + ω_{i + 1}^{2} r_{i + 1} \end{array}

\begin{array}{ll} \frac{d^{2}}{d t^{2}} (d_{i} + r_{i}) & = \frac{d^{2}}{d t^{2}} (d_{i + 1} - r_{i + 1}) \\ = a_{i} + α_{i} \times r_{i} - ω_{i}^{2} r_{i} \\ = a_{i + 1} - α_{i + 1} \times r_{i + 1} + ω_{i + 1}^{2} r_{i + 1} \end{array}

Multiplying both sides by $m$ gives:

m (a_{i} + α_{i} \times r_{i}) - m ω_{i}^{2} r_{i} = m (a_{i + 1} - α_{i + 1} \times r_{i + 1}) + m ω_{i + 1}^{2} r_{i + 1}

If $F_{res, i}$ is the total force applied to i-rod other than $F_{i}$ and $- F_{i - 1}$ , which is $F_{res, i} = F_{c κ, i} + F_{b, i}$ , then:

m a_{i} = F_{i} - F_{i - 1} + F_{res, i}

If $τ_{res, i}$ is the total torque applied to i-rod excluding $r_{i} \times (F_{i - 1} + F_{i})$ , which is $τ_{res, i} = τ_{c κ, i} + τ_{b, i}$ , then:

I α_{i} = \frac{1}{3} m r^{2} α_{i} = r_{i} \times (F_{i - 1} + F_{i}) + τ_{res, i}

Dividing both sides by $\frac{1}{3} r^{2}$ gives:

m α_{i} = 3 r^{- 2} r_{i} \times (F_{i - 1} + F_{i}) + 3 r^{- 2} τ_{res, i}

If $h_{res, i} \equiv 3 r^{- 2} τ_{res, i} \times r_{i}$ , then:

\begin{array}{ll} m (a_{i} + α_{i} \times r_{i}) & = (F_{i} - F_{i - 1} + F_{res, i}) + [3 r^{- 2} r_{i} \times (F_{i - 1} + F_{i}) + 3 r^{- 2} τ_{res, i}] \times r_{i} \\ = (F_{i} - F_{i - 1}) + 3 [{\hat{r}}_{i} \times (F_{i - 1} + F_{i})] \times {\hat{r}}_{i} + F_{res, i} + 3 r^{- 2} τ_{res, i} \times r_{i} \\ = (F_{i} - F_{i - 1}) + 3 [(F_{i - 1} + F_{i}) \cdot {\hat{N}}_{i}] {\hat{N}}_{i} + F_{res, i} + h_{res, i} \end{array}

Following the same method,

m (a_{i + 1} - α_{i + 1} \times r_{i + 1}) = (F_{i + 1} - F_{i}) - 3 [(F_{i} + F_{i + 1}) \cdot {\hat{N}}_{i + 1}] {\hat{N}}_{i + 1} + F_{res, i + 1} - h_{res, i + 1}

To organize the terms of the equations into known values and unknown values,

\begin{aligned} k_{i} & \equiv F_{i} - F_{i - 1} \\ P_{i} & \equiv 3 {\hat{N}}_{i} {\hat{N}}_{i}^{T} \\ h_{i} & \equiv 3 [(F_{i - 1} + F_{i}) \cdot {\hat{N}}_{i}] {\hat{N}}_{i} \\ = 3 {\hat{N}}_{i} {\hat{N}}_{i}^{T} (F_{i - 1} + F_{i}) \\ = P_{i} (F_{i - 1} + F_{i}) \\ m (a_{i} + α_{i} \times r_{i}) & = k_{i} + h_{i} + F_{res, i} + h_{res, i} \\ m (a_{i + 1} - α_{i + 1} \times r_{i + 1}) & = k_{i + 1} - h_{i + 1} + F_{res, i + 1} - h_{res, i + 1} \\ k_{i} + h_{i} + F_{res, i} + h_{res, i} - m ω_{i}^{2} r_{i} & = k_{i + 1} - h_{i + 1} + F_{res, i + 1} - h_{res, i + 1} + m ω_{i + 1}^{2} r_{i + 1} \\ q_{i} & \equiv k_{i} + h_{i} - k_{i + 1} + h_{i + 1} \\ = - F_{res, i} - h_{res, i} + F_{res, i + 1} - h_{res, i + 1} + m ω_{i}^{2} r_{i} + m ω_{i + 1}^{2} r_{i + 1} \end{aligned}

Therefore, if we know all $F_{res, i}$ and $τ_{res, i}$ for each i, we can find $q_{i}$ . If $I \equiv [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}]$ , and if we expand $q_{i}$ ,

\begin{array}{ll} q_{i} & \equiv k_{i} + h_{i} - k_{i + 1} + h_{i + 1} \\ = - k_{i + 1} + k_{i} + h_{i} + h_{i + 1} \\ = - F_{i + 1} + F_{i} + F_{i} - F_{i - 1} + P_{i} (F_{i - 1} + F_{i}) + P_{i + 1} (F_{i} + F_{i + 1}) \\ = F_{i - 1} (P_{i} - I) + F_{i} (P_{i} + P_{i + 1} + 2 I) + F_{i + 1} (P_{i + 1} - I) \end{array}

If we set $A_{i} \equiv P_{i} - I$ and $B_{i} \equiv P_{i} + P_{i + 1} + 2 I$ , then:

q_{i} = A_{i} F_{i - 1} + B_{i} F_{i} + A_{i + 1} F_{i + 1}

If we set $\bar{0} = [\begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix}]$ and express the above equation in a multidimensional tensor form,

[\begin{matrix} B_{1} & A_{2} & \bar{0} & \bar{0} & \bar{0} & \dots & \bar{0} \\ A_{2} & B_{2} & A_{3} & \bar{0} & \bar{0} & \dots & \bar{0} \\ \bar{0} & A_{3} & B_{3} & A_{4} & \bar{0} & \dots & \bar{0} \\ \bar{0} & \bar{0} & A_{4} & B_{4} & A_{5} & \dots & \bar{0} \\ \bar{0} & \bar{0} & \bar{0} & A_{5} & B_{5} & ⋱ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & A_{n - 1} \\ \bar{0} & \bar{0} & \bar{0} & \bar{0} & \dots & A_{n - 1} & B_{n - 1} \end{matrix}] [\begin{matrix} F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ ⋮ \\ F_{n - 2} \\ F_{n - 1} \end{matrix}] = [\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ ⋮ \\ q_{n - 2} \\ q_{n - 1} \end{matrix}]

Let us set:

D \equiv [\begin{matrix} B_{1} & A_{2} & \bar{0} & \bar{0} & \bar{0} & \dots & \bar{0} \\ A_{2} & B_{2} & A_{3} & \bar{0} & \bar{0} & \dots & \bar{0} \\ \bar{0} & A_{3} & B_{3} & A_{4} & \bar{0} & \dots & \bar{0} \\ \bar{0} & \bar{0} & A_{4} & B_{4} & A_{5} & \dots & \bar{0} \\ \bar{0} & \bar{0} & \bar{0} & A_{5} & B_{5} & ⋱ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & A_{n - 1} \\ \bar{0} & \bar{0} & \bar{0} & \bar{0} & \dots & A_{n - 1} & B_{n - 1} \end{matrix}]

F \equiv [\begin{matrix} F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ ⋮ \\ F_{n - 2} \\ F_{n - 1} \end{matrix}]

Q \equiv [\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ ⋮ \\ q_{n - 2} \\ q_{n - 1} \end{matrix}]

As $A_{i}$ and $B_{i}$ can be calculated from $P_{i}$ , $P_{i}$ from ${\hat{N}}_{i}$ , and ${\hat{N}}_{i}$ from $s_{i}$ , we can find $D$ from $s_{i}$ .

\begin{array}{ll} q_{i} & = - F_{res, i} + F_{res, i + 1} - h_{res, i} - h_{res, i + 1} + m ω_{i}^{2} r_{i} + m ω_{i + 1}^{2} r_{i + 1} \\ = - F_{res, i} + F_{res, i + 1} - 3 r^{- 2} (τ_{res, i} \times r_{i} + τ_{res, i + 1} \times r_{i + 1}) + m ω_{i}^{2} r_{i} + m ω_{i + 1}^{2} r_{i + 1} \end{array}

Therefore, if we know $F_{c κ, i}$ , $F_{b, i}$ , $τ_{c κ, i}$ , $τ_{b, i}$ , we can find $Q$ .

∴ D F = Q

We can find $F$ , thus $F_{i}$ ( $i \in {1, \dots, n - 1}$ ), by solving this tensor equation.

If we denote each component of the tensor, which are each matrix and the p-th row q-th column component of the vector, as $(*)_{p, q}$ , then:

[\begin{matrix} [\begin{matrix} (B_{1})_{1, 1} & (B_{1})_{1, 2} \\ (B_{1})_{2, 1} & (B_{1})_{2, 2} \end{matrix}] & [\begin{matrix} (A_{2})_{1, 1} & (A_{2})_{1, 2} \\ (A_{2})_{2, 1} & (A_{2})_{2, 2} \end{matrix}] & [\begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix}] & \dots \\ [\begin{matrix} (A_{2})_{1, 1} & (A_{2})_{1, 2} \\ (A_{2})_{2, 1} & (A_{2})_{2, 2} \end{matrix}] & [\begin{matrix} (B_{2})_{1, 1} & (B_{2})_{1, 2} \\ (B_{2})_{2, 1} & (B_{2})_{2, 2} \end{matrix}] & [\begin{matrix} (A_{3})_{1, 1} & (A_{3})_{1, 2} \\ (A_{3})_{2, 1} & (A_{3})_{2, 2} \end{matrix}] & \dots \\ [\begin{matrix} 0 & 0 \\ 0 & 0 \end{matrix}] & [\begin{matrix} (A_{3})_{1, 1} & (A_{3})_{1, 2} \\ (A_{3})_{2, 1} & (A_{3})_{2, 2} \end{matrix}] & [\begin{matrix} (B_{3})_{1, 1} & (B_{3})_{1, 2} \\ (B_{3})_{2, 1} & (B_{3})_{2, 2} \end{matrix}] & \dots \\ ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}] [\begin{matrix} [\begin{matrix} (F_{1})_{1, 1} \\ (F_{1})_{2, 1} \end{matrix}] \\ [\begin{matrix} (F_{2})_{1, 1} \\ (F_{2})_{2, 1} \end{matrix}] \\ [\begin{matrix} (F_{3})_{1, 1} \\ (F_{3})_{2, 1} \end{matrix}] \\ ⋮ \end{matrix}] = [\begin{matrix} [\begin{matrix} (q_{1})_{1, 1} \\ (q_{1})_{2, 1} \end{matrix}] \\ [\begin{matrix} (q_{2})_{1, 1} \\ (q_{2})_{2, 1} \end{matrix}] \\ [\begin{matrix} (q_{3})_{1, 1} \\ (q_{3})_{2, 1} \end{matrix}] \\ ⋮ \end{matrix}]

And the matrix equation equivalent to this tensor equation is:

[\begin{matrix} (B_{1})_{1, 1} & (B_{1})_{1, 2} & (A_{2})_{1, 1} & (A_{2})_{1, 2} & 0 & 0 & \dots \\ (B_{1})_{2, 1} & (B_{1})_{2, 2} & (A_{2})_{2, 1} & (A_{2})_{2, 2} & 0 & 0 & \dots \\ (A_{2})_{1, 1} & (A_{2})_{1, 2} & (B_{2})_{1, 1} & (B_{2})_{1, 2} & (A_{3})_{1, 1} & (A_{3})_{1, 2} & \dots \\ (A_{2})_{2, 1} & (A_{2})_{2, 2} & (B_{2})_{2, 1} & (B_{2})_{2, 2} & (A_{3})_{2, 1} & (A_{3})_{2, 2} & \dots \\ 0 & 0 & (A_{3})_{1, 1} & (A_{3})_{1, 2} & (B_{3})_{1, 1} & (B_{3})_{1, 2} & \dots \\ 0 & 0 & (A_{3})_{2, 1} & (A_{3})_{2, 2} & (B_{3})_{2, 1} & (B_{3})_{2, 2} & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}] [\begin{matrix} (F_{1})_{1, 1} \\ (F_{1})_{2, 1} \\ (F_{2})_{1, 1} \\ (F_{2})_{2, 1} \\ (F_{3})_{1, 1} \\ (F_{3})_{2, 1} \\ ⋮ \end{matrix}] = [\begin{matrix} (q_{1})_{1, 1} \\ (q_{1})_{2, 1} \\ (q_{2})_{1, 1} \\ (q_{2})_{2, 1} \\ (q_{3})_{1, 1} \\ (q_{3})_{2, 1} \\ ⋮ \end{matrix}]

Let us set:

D \equiv [\begin{matrix} (B_{1})_{1, 1} & (B_{1})_{1, 2} & (A_{2})_{1, 1} & (A_{2})_{1, 2} & 0 & 0 & \dots \\ (B_{1})_{2, 1} & (B_{1})_{2, 2} & (A_{2})_{2, 1} & (A_{2})_{2, 2} & 0 & 0 & \dots \\ (A_{2})_{1, 1} & (A_{2})_{1, 2} & (B_{2})_{1, 1} & (B_{2})_{1, 2} & (A_{3})_{1, 1} & (A_{3})_{1, 2} & \dots \\ (A_{2})_{2, 1} & (A_{2})_{2, 2} & (B_{2})_{2, 1} & (B_{2})_{2, 2} & (A_{3})_{2, 1} & (A_{3})_{2, 2} & \dots \\ 0 & 0 & (A_{3})_{1, 1} & (A_{3})_{1, 2} & (B_{3})_{1, 1} & (B_{3})_{1, 2} & \dots \\ 0 & 0 & (A_{3})_{2, 1} & (A_{3})_{2, 2} & (B_{3})_{2, 1} & (B_{3})_{2, 2} & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}]

Since $A_{i}$ , $B_{i}$ are symmetric matrices, $D$ is a heptadiagonal symmetric matrix. Since $D$ is a symmetric matrix, the solution to the matrix equation can be found with the Cholesky decomposition. Therefore, we can find the x-axis and y-axis components of $F_{i}$ .

As a result, we can find the joint force $F_{i}$ from the known values ( $s_{i}$ , $F_{c κ, i}$ , $F_{b, i}$ , $τ_{c κ, i}$ , $τ_{b, i}$ ).

Proof of numerical integration for the translational motion of a worm using semi-implicit Euler method

When the friction coefficients $b_{⊥}$ , $b_{∥}$ are sufficiently large compared to $M / Δ t$ , numerical integration via the explicit Euler method ( $v_{c}^{(t + Δ t)} = v_{c}^{(t)} + a_{c}^{(t)} Δ t = v_{c}^{(t)} + \frac{\sum_{i} F_{b, i}^{(t)}}{M} Δ t$ ) becomes unstable (Butcher, 2003). Therefore, for all friction coefficients greater than or equal to 0, the semi-implicit Euler method ( $v_{c}^{(t + Δ t)} = v_{c}^{(t)} + a_{c}^{(t + Δ t)} Δ t ≃ v_{c}^{(t)} + \frac{1}{1 + \frac{b_{⊥} Δ t}{M}} \frac{\sum_{i} F_{b, i}^{(t)}}{M} Δ t$ ) was used to ensure numerical integration remains stable, and its proof is as follows.

Newton’s equation for the translational motion of each i-rod is as follows.

\begin{array}{ll} F_{i} & = m a_{i} \\ = F_{b, i} + F_{c κ, i} + F_{j o i n t, i} \\ = - \frac{b_{⊥}}{n} v_{⊥, i} - \frac{b_{∥}}{n} v_{∥, i} + F_{c κ, i} + F_{j o i n t, i} \end{array}

The integration formula for $a_{i}$ using the implicit Euler method is as follows. (where $b_{⊥} \geq b_{∥}$ )

\begin{aligned} v_{i}^{(t + Δ t)} & = v_{i}^{(t)} + a_{i}^{(t + Δ t)} Δ t \\ = v_{i}^{(t)} - \frac{Δ t}{m} (\frac{b_{⊥}}{n} v_{⊥, i}^{(t + Δ t)} + \frac{b_{∥}}{n} v_{∥, i}^{(t + Δ t)}) + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ = v_{i}^{(t)} - \frac{Δ t}{M} (b_{⊥} v_{⊥, i}^{(t + Δ t)} + b_{∥} v_{∥, i}^{(t + Δ t)}) + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ = v_{i}^{(t)} - \frac{Δ t}{M} (b_{⊥} (v_{⊥, i}^{(t + Δ t)} + v_{∥, i}^{(t + Δ t)}) + (b_{∥} - b_{⊥}) v_{∥, i}^{(t + Δ t)}) + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ = v_{i}^{(t)} - \frac{Δ t}{M} (b_{⊥} v_{i}^{(t + Δ t)} + (b_{∥} - b_{⊥}) v_{∥, i}^{(t + Δ t)}) + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ (1 + \frac{b_{⊥} Δ t}{M}) v_{i}^{(t + Δ t)} & = v_{i}^{(t)} - \frac{Δ t}{M} ((b_{∥} - b_{⊥}) v_{∥, i}^{(t + Δ t)}) + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ = v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t + Δ t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \end{aligned}

Because $v_{⊥, i}^{(t + Δ t)}$ and $v_{∥, i}^{(t + Δ t)}$ are unknown at time $t$ , the numerical calculation of the above formula is impossible.

\frac{| \frac{(b_{⊥} - b_{∥}) Δ t}{M} (v_{∥, i}^{(t + Δ t)} - v_{∥, i}^{(t)}) |}{| v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t + Δ t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) |} ≃ 0

If the above formula is assumed to be true, the following approximation can be used.

\begin{aligned} (1 + \frac{b_{⊥} Δ t}{M}) v_{i}^{(t + Δ t)} & = v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t + Δ t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \\ ≃ v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}) \end{aligned}

The approximate value of $v_{i}^{(t + Δ t)}$ by the above approximation is as follows.

v_{i}^{(t + Δ t)} ≃ [v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)})] / (1 + \frac{b_{⊥} Δ t}{M})

The approximate value of $a_{i}^{(t + Δ t)} Δ t$ by the approximate value of $v_{i}^{(t + Δ t)}$ is as follows.

\begin{array}{ll} a_{i}^{(t + Δ t)} Δ t & = v_{i}^{(t + Δ t)} - v_{i}^{(t)} \\ ≃ [v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)})] / (1 + \frac{b_{⊥} Δ t}{M}) - v_{i}^{(t)} \\ = [- \frac{b_{⊥} Δ t}{M} v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥}) Δ t}{M} v_{∥, i}^{(t)} + \frac{Δ t}{m} (F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)})] / (1 + \frac{b_{⊥} Δ t}{M}) \\ = \frac{Δ t}{m} [- \frac{b_{⊥}}{n} v_{i}^{(t)} + \frac{(b_{⊥} - b_{∥})}{n} v_{∥, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}] / (1 + \frac{b_{⊥} Δ t}{M}) \\ = \frac{Δ t}{m} [- \frac{b_{⊥}}{n} (v_{i}^{(t)} - v_{∥, i}^{(t)}) - \frac{b_{∥}}{n} v_{∥, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}] / (1 + \frac{b_{⊥} Δ t}{M}) \\ = \frac{Δ t}{m} [- \frac{b_{⊥}}{n} v_{⊥, i}^{(t)} - \frac{b_{∥}}{n} v_{∥, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}] / (1 + \frac{b_{⊥} Δ t}{M}) \\ = \frac{Δ t}{m} [F_{b, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}] / (1 + \frac{b_{⊥} Δ t}{M}) \end{array}

If both sides are divided by $Δ t$ ,

a_{i}^{(t + Δ t)} ≃ \frac{1}{(1 + \frac{b_{⊥} Δ t}{M})} \frac{F_{b, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}}{m}

$F_{c κ, i}^{(t + Δ t)}$ , $F_{j o i n t, i}^{(t + Δ t)}$ satisfy the following because they are the internal forces of the worm ( $\vec{0}$ is a zero vector).

\begin{aligned} \sum_{i} F_{c κ, i}^{(t + Δ t)} & = \vec{0} \\ \sum_{i} F_{j o i n t, i}^{(t + Δ t)} & = \vec{0} \end{aligned}

Therefore, the approximation of the force $M a_{c}^{(t + Δ t)}$ received by the worm is as follows.

\begin{array}{ll} M a_{c}^{(t + Δ t)} & = \sum_{i} m a_{i}^{(t + Δ t)} \\ ≃ \sum_{i} m \frac{1}{(1 + \frac{b_{⊥} Δ t}{M})} \frac{F_{b, i}^{(t)} + F_{c κ, i}^{(t + Δ t)} + F_{j o i n t, i}^{(t + Δ t)}}{m} \\ = \frac{1}{(1 + \frac{b_{⊥} Δ t}{M})} (\sum_{i} F_{b, i}^{(t)} + \sum_{i} F_{c κ, i}^{(t + Δ t)} + \sum_{i} F_{j o i n t, i}^{(t + Δ t)}) \\ = \frac{1}{(1 + \frac{b_{⊥} Δ t}{M})} \sum_{i} F_{b, i}^{(t)} \end{array}

If both sides are divided by $M$ ,

a_{c}^{(t + Δ t)} ≃ \frac{1}{(1 + \frac{b_{⊥} Δ t}{M})} \frac{\sum_{i} F_{b, i}^{(t)}}{M}

This approximation ensures computational stability regardless of the size of $b_{⊥}$ , $b_{∥}$ . That is, this approximation solves the problem of the decrease in computational stability of numerical integration through the explicit Euler method when $b_{⊥}$ , $b_{∥}$ are sufficiently large compared to $M / Δ t$ .

Numerical integration of the rotational motion of i-rod using semi-implicit Euler method

First, the numerical integration formulas for $ω_{i}$ and $α_{i}$ using the implicit Euler method are as follows.

s_{i}^{(t + Δ t)} = s_{i}^{(t)} + ω_{i}^{(t + Δ t)} Δ t

ω_{i}^{(t + Δ t)} = ω_{i}^{(t)} + α_{i}^{(t + Δ t)} Δ t

If we set $β \equiv \frac{1}{3} \frac{b_{⊥}}{n} r^{2}$ , the equation describing the rotation of i-rod is as follows.

\begin{aligned} I α_{i} & = τ_{b, i} + τ_{c κ, i} + τ_{j o i n t, i} \\ = - β ω_{i} + c (ω_{i + 1} - 2 ω_{i} + ω_{i - 1}) + κ (θ_{i} - θ_{i - 1} - θ_{c t r l, i} + θ_{c t r l, i - 1}) + τ_{j o i n t, i} \\ = - β ω_{i} + c (ω_{i + 1} - 2 ω_{i} + ω_{i - 1}) + κ (s_{i + 1} - 2 s_{i} + s_{i - 1} - θ_{c t r l, i} + θ_{c t r l, i - 1}) + τ_{j o i n t, i} \\ = c ω_{i - 1} - (β + 2 c) ω_{i} + c ω_{i + 1} + κ (s_{i - 1} - 2 s_{i} + s_{i + 1}) - κ (θ_{c t r l, i} - θ_{c t r l, i + 1}) + τ_{j o i n t, i} \end{aligned}

If the above formula is expanded for time $t + Δ t$ ,

\begin{aligned} I α_{i}^{(t + Δ t)} = & c ω_{i - 1}^{(t + Δ t)} - (β + 2 c) ω_{i}^{(t + Δ t)} + c ω_{i + 1}^{(t + Δ t)} \\ + κ (s_{i - 1}^{(t + Δ t)} - 2 s_{i}^{(t + Δ t)} + s_{i + 1}^{(t + Δ t)}) - κ (θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t + Δ t)}) + τ_{j o i n t, i}^{(t + Δ t)} \\ = & c ω_{i - 1}^{(t + Δ t)} - (β + 2 c) ω_{i}^{(t + Δ t)} + c ω_{i + 1}^{(t + Δ t)} \\ + κ ((s_{i - 1}^{(t)} + ω_{i - 1}^{(t + Δ t)} Δ t) - 2 (s_{i}^{(t)} + ω_{i}^{(t + Δ t)} Δ t) + (s_{i + 1}^{(t)} + ω_{i + 1}^{(t + Δ t)} Δ t)) \\ - κ (θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t + Δ t)}) + τ_{j o i n t, i}^{(t + Δ t)} \\ = & (c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} \\ + κ (s_{i - 1}^{(t)} - 2 s_{i}^{(t)} + s_{i + 1}^{(t)}) - κ (θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t + Δ t)}) + τ_{j o i n t, i}^{(t + Δ t)} \end{aligned}

The above formula is impossible to integrate because $θ_{c t r l, i}^{(t + Δ t)}$ , $θ_{c t r l, i + 1}^{(t + Δ t)}$ , $τ_{j o i n t, i}^{(t + Δ t)}$ are unknown at time t.

\frac{| - κ ((θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i}^{(t)}) - (θ_{c t r l, i + 1}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t)})) + (τ_{j o i n t, i}^{(t + Δ t)} - τ_{j o i n t, i}^{(t)}) |}{| (c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} + κ (s_{i - 1}^{(t)} - 2 s_{i}^{(t)} + s_{i + 1}^{(t)}) - κ (θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t + Δ t)}) + τ_{j o i n t, i}^{(t + Δ t)} |} ≃ 0

If the above formula is assumed to be true, the following approximation can be used.

\begin{aligned} I α_{i}^{(t + Δ t)} = & (c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} \\ + κ (s_{i - 1}^{(t)} - 2 s_{i}^{(t)} + s_{i + 1}^{(t)}) - κ (θ_{c t r l, i}^{(t + Δ t)} - θ_{c t r l, i + 1}^{(t + Δ t)}) + τ_{j o i n t, i}^{(t + Δ t)} \\ ≃ & (c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} \\ + κ (s_{i - 1}^{(t)} - 2 s_{i}^{(t)} + s_{i + 1}^{(t)}) - κ (θ_{c t r l, i}^{(t)} - θ_{c t r l, i + 1}^{(t)}) + τ_{j o i n t, i}^{(t)} \\ = & (c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} + τ_{κ, i}^{(t)} + τ_{j o i n t, i}^{(t)} \\ ω_{i}^{(t + Δ t)} = & ω_{i}^{(t)} + α_{i}^{(t + Δ t)} Δ t \\ = & ω_{i}^{(t)} + \frac{Δ t}{I} ((c + κ Δ t) ω_{i - 1}^{(t + Δ t)} - (β + 2 c + 2 κ Δ t) ω_{i}^{(t + Δ t)} + (c + κ Δ t) ω_{i + 1}^{(t + Δ t)} + τ_{κ, i}^{(t)} + τ_{j o i n t, i}^{(t)}) \end{aligned}

The above formula can be expressed as a matrix formula as follows.

[\begin{matrix} ⋮ \\ ω_{i - 1}^{(t + Δ t)} \\ ω_{i}^{(t + Δ t)} \\ ω_{i + 1}^{(t + Δ t)} \\ ⋮ \end{matrix}] ≃ [\begin{matrix} ⋮ \\ ω_{i - 1}^{(t)} \\ ω_{i}^{(t)} \\ ω_{i + 1}^{(t)} \\ ⋮ \end{matrix}] + \frac{Δ t}{I} [\begin{matrix} ⋱ & ⋮ & ⋮ & ⋮ & ⋱ \\ \dots & (c + κ Δ t) & - (β + 2 c + 2 κ Δ t) & (c + κ Δ t) & \dots \\ ⋱ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}] [\begin{matrix} ⋮ \\ ω_{i - 1}^{(t + Δ t)} \\ ω_{i}^{(t + Δ t)} \\ ω_{i + 1}^{(t + Δ t)} \\ ⋮ \end{matrix}] + \frac{Δ t}{I} [\begin{matrix} ⋮ \\ τ_{κ, i - 1}^{(t)} + τ_{j o i n t, i - 1}^{(t)} \\ τ_{κ, i}^{(t)} + τ_{j o i n t, i}^{(t)} \\ τ_{κ, i + 1}^{(t)} + τ_{j o i n t, i + 1}^{(t)} \\ ⋮ \end{matrix}]

In the above vector matrix formula, let us represent the vectors and matrix by the following symbols.

{\vec{ω}}^{(t)} \equiv [\begin{matrix} ⋮ \\ ω_{i - 1}^{(t)} \\ ω_{i}^{(t)} \\ ω_{i + 1}^{(t)} \\ ⋮ \end{matrix}]

P_{n \times n} \equiv [\begin{matrix} ⋱ & ⋮ & ⋮ & ⋮ & ⋱ \\ \dots & (c + κ Δ t) & - (β + 2 c + 2 κ Δ t) & (c + κ Δ t) & \dots \\ ⋱ & ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}]

{\vec{τ}}_{r e m}^{(t)} \equiv [\begin{matrix} ⋮ \\ τ_{κ, i - 1}^{(t)} + τ_{j o i n t, i - 1}^{(t)} \\ τ_{κ, i}^{(t)} + τ_{j o i n t, i}^{(t)} \\ τ_{κ, i + 1}^{(t)} + τ_{j o i n t, i + 1}^{(t)} \\ ⋮ \end{matrix}]

Then, the matrix formula is expressed as follows.

{\vec{ω}}^{(t + Δ t)} ≃ {\vec{ω}}^{(t)} + \frac{Δ t}{I} P_{n \times n} {\vec{ω}}^{(t + Δ t)} + \frac{Δ t}{I} {\vec{τ}}_{r e m}^{(t)}

Now, the following approximation can be obtained where $I_{n \times n}$ is a unit matrix of size $n \times n$ .

∴ {\vec{ω}}^{(t + Δ t)} ≃ {(I_{n \times n} - \frac{Δ t}{I} P_{n \times n})}^{- 1} ({\vec{ω}}^{(t)} + \frac{Δ t}{I} {\vec{τ}}_{r e m}^{(t)})

Correction formula for the rotational inertia of the entire worm

For a floor surface with low friction like water, when numerically integrating the rotational motion of the worm, if $Δ t > 1 \times 10^{- 6} s e c$ , the calculation error accumulated for the rotational inertia of the whole worm significantly influenced the calculation result $ω_{i}$ and $s_{i}$ . To prevent this, the error is corrected as follows. If ${\bar{x}}_{i} \equiv x_{i} - x_{c}$ , ${\bar{y}}_{i} \equiv y_{i} - y_{c}$ , the moment of inertia of the entire worm at time t is as follows by the parallel axis theorem.

I_{b o d y}^{(t)} = m \sum_{i} ({({\bar{x}}_{i}^{(t)})}^{2} + {({\bar{y}}_{i}^{(t)})}^{2}) + n I

If i-rod is approximated as a point particle, the torque applied to the entire worm at time t is as follows where subscription x, y indicates x, y components of the vector. (See ‘Numerical integration for translational motion’ in Appendix)

τ_{b o d y}^{(t + Δ t)} ≃ \frac{\sum_{i} ({\bar{x}}_{i}^{(t)} F_{b, i, y}^{(t)} - {\bar{y}}_{i}^{(t)} F_{b, i, x}^{(t)})}{1 + \frac{b_{⊥} Δ t}{M}}

If i-rod is approximated as a point particle, the rotational inertia of the whole worm at time $t$ is as follows.

L_{b o d y}^{(t)} ≃ m \sum_{i} ({\bar{x}}_{i}^{(t)} v_{i, y}^{(t)} - {\bar{y}}_{i}^{(t)} v_{i, x}^{(t)})

The predicted value of $ω_{i}$ at $t + Δ t$ , $ω_{i}^{p}$ , is calculated by the semi-implicit Euler method (See ‘Numerical integration of the rotational motion’ in Appendix). ( ${\vec{ω}}^{p} = {[\begin{matrix} \dots & ω_{i}^{p} & \dots \end{matrix}]}^{T}$ )

{\vec{ω}}^{p} ≃ {(I_{n \times n} - \frac{Δ t}{I} P_{n \times n})}^{- 1} ({\vec{ω}}^{(t)} + \frac{Δ t}{I} {\vec{τ}}_{r e m}^{(t)})

The predicted value of $s_{i}$ at $t + Δ t$ is as follows.

s_{i}^{p} = s_{i}^{(t)} + ω_{i}^{p} Δ t

Predicted values $x_{i}^{p}$ , $y_{i}^{p}$ , $v_{i}^{p}$ for $x_{i}^{(t + Δ t)}$ , $y_{i}^{(t + Δ t)}$ , $v_{i}^{(t + Δ t)}$ are calculated from $d_{c}$ , $v_{c}$ , $s_{i}^{p}$ , $ω_{i}^{p}$ (See ‘Minimum information required to describe the motion of each rod’ in Appendix).

Using $x_{i}^{p}$ , $y_{i}^{p}$ , $v_{i}^{p}$ , the moment of inertia $I_{b o d y}^{p}$ and rotational inertia $L_{b o d y}^{p}$ of the entire worm at time $t + Δ t$ are calculated as follows. (where ${\bar{x}}_{i}^{p} \equiv x_{i}^{p} - x_{c}^{p}$ , ${\bar{y}}_{i}^{p} \equiv y_{i}^{p} - y_{c}^{p}$ , $v_{i} = {[\begin{matrix} v_{i, x} & v_{i, y} \end{matrix}]}^{T}$ )

I_{b o d y}^{p} = m \sum_{i} ({({\bar{x}}_{i}^{p})}^{2} + {({\bar{y}}_{i}^{p})}^{2}) + n I

L_{b o d y}^{p} = m \sum_{i} ({\bar{x}}_{i}^{p} v_{i, y}^{p} - {\bar{y}}_{i}^{p} v_{i, x}^{p})

$ω_{i}^{(t + Δ t)}$ is calculated as follows.

ω_{i}^{(t + Δ t)} = ω_{i}^{p} + \frac{- (L_{b o d y}^{p} - L_{b o d y}^{(t)}) + τ_{b o d y}^{(t + Δ t)} Δ t}{\frac{I_{b o d y}^{p} + I_{b o d y}^{(t)}}{2}}

$s_{i}^{(t + Δ t)}$ is calculated as follows.

s_{i}^{(t + Δ t)} = s_{i}^{(t)} + ω_{i}^{(t + Δ t)} Δ t

By correcting the rotational inertia for the whole worm, numerical integration of the rotational motion of the worm was well calculated even for cases when $Δ t > 1 \times 10^{- 6} s e c$ , as if $Δ t \leq 1 \times 10^{- 6} s e c$ .

Proper selection of friction coefficients

When using the vertical and horizontal friction coefficients $b_{a g a r, ⊥}$ , $b_{a g a r, ∥}$ on agar, as proposed in the previous work (Boyle et al., 2012), the trajectory of the escaping behavior was not accurately replicated. Therefore, we sought appropriate friction coefficients necessary for replicating the escaping behavior. We calculated new vertical and horizontal friction coefficients $b_{η, ⊥} = η b_{a g a r, ⊥}$ , $b_{η, ∥} = η b_{a g a r, ∥}$ by multiplying scaling factor $η$ to the agar friction coefficients $b_{a g a r, ⊥}$ , $b_{a g a r, ∥}$ of the previous work (Boyle et al., 2012). Let us denote the set of a quantity for all pairs $(i, t)$ of index i and time $t$ as ${*}_{i, t}$ . When the kymogram input ${θ_{c t r l, i}^{(t)}}_{i, t}$ was same as Figure 3A, we observed how the trajectory of the escaping behavior changes with $η$ (Appendix 1—figure 3A) and analyzed representative values for each trajectory as follows (Appendix 1—figure 3B) (Note that when $η$ was less than or equal to 10^-6 , the time-step ( $Δ t$ ) was set to 10^-6 sec for higher accuracy of simulation). Let us denote the average of a quantity for all pairs $(i, t)$ as $⟨ * ⟩_{i, t}$ . When $η$ was 1 ∼ 10^-9 , the smaller $η$ , the smaller $E_{θ} = {⟨ | θ_{i}^{(t)} - θ_{c t r l, i}^{(t)} | ⟩}_{i, t}$ was. The reduction in $E_{θ}$ when $η$ changed from 10^-2 to 10^-9 was about 3% of the reduction when $η$ changed from 1 to 10^-9 . When $η$ was between 1 and 10^-9 , even if $η$ decreased, $E_{θ}$ did not fall below 0.125 (rad). This is because there was a time delay between the input $θ_{c t r l, i}^{(t)}$ and the response $θ_{i}^{(t)}$ . As $η$ decreased from 1 to 10^-2 , the total traveled distance of the worm ( $\sum_{t} | v_{c}^{(t)} Δ t |$ ) and the total absolute angle change ( $S = \sum_{t = 0}^{T - Δ t} | {⟨ s_{i}^{(t + Δ t)} ⟩}_{i} - {⟨ s_{i}^{(t)} ⟩}_{i} |$ ) increased, and the trajectory became more similar to the experimental video (Broekmans et al., 2016). When $η$ was between 10^-2 and 10^-6 , the worm’s trajectory was almost identical, and thus the total traveled distance, $S$ , and the pattern of $| F_{b, i}^{(t)} |$ (Appendix 1—figure 4) were similar across trajectories. If $η$ was smaller than 10^-6 , the worm’s total traveled distance decreased and $S$ increased, and the trajectory was no longer similar to the experimental video. This was because a too small friction coefficients hindered the worm from obtaining enough propulsive force from the ground (Appendix 1—figure 4).

The same method with the kymogram input ${θ_{c t r l, i}^{(t)}}_{i, t}$ same as Figure 3B used to analyze the effect of $η$ on the trajectory of the escaping behavior was applied to analyze the trajectory of the delta-turn according to $η$ (Appendix 1—figure 5A). When $η$ was 1 ∼ 10^-9 , the smaller $η$ , the smaller $E_{θ}$ was (Appendix 1—figure 5B). The reduction in $E_{θ}$ when $η$ changed from 10^-2 to 10^-9 was about 3% of the reduction when $η$ changed from 1 to 10^-9 . When $η$ was between 1 and 10^-9 , even if $η$ decreased, $E_{θ}$ did not fall below 0.15 (rad). As $η$ decreased from 1 to 10^-2 , $S$ increased, and the trajectory became more similar to the experimental video. When $η$ was between 10^-2 and 10^-6 , the worm’s trajectory was almost identical, and thus the total traveled distance, $S$ , and the pattern of $| F_{b, i}^{(t)} |$ (Appendix 1—figure 6) were similar across trajectories. If $η$ was smaller than 10^-6 , the worm’s total traveled distance decreased and $S$ increased, and the trajectory was no longer similar to the experimental video. This is because a too small friction coefficients hindered the worm from obtaining enough propulsive force from the ground (Appendix 1—figure 6).

In conclusion, when the ratio between vertical and horizontal friction coefficients was constant at 40( $= b_{η, ⊥} / b_{η, ∥} = b_{a g a r, ⊥} / b_{a g a r, ∥}$ ), selecting an appropriate $η$ value (between 10^-2 and 10^-6) was crucial for replicating the trajectory of sequenced locomotive behavior. Therefore, we chose $b_{η, ⊥}$ , $b_{η, ∥}$ ( $η = 10^{- 2}$ ) as the friction coefficients that sufficiently reduce $E_{θ}$ among those closest to the agar friction coefficients of the previous work (Boyle et al., 2012) for both escaping behavior and delta-turn.

Appendix 1—figure 1

Download asset Open asset

Method of compressing motion state information.

(A) Method of calculating the relative position of the rod. (B) Method of calculating the absolute position of the rod. (C) Method of calculating the relative velocity of the rod.

Appendix 1—figure 2

Download asset Open asset

The i-actuator.

(A) Composition of i-actuator. (B) Forces that i-actuator gives to i-rod and (i+1)-rod. (C) Resultant force given by i-actuator. (D) Force component of i-actuator that applies torque to i-rod.

Appendix 1—figure 3

Download asset Open asset

The effect of the scaling factor η of the friction coefficients on the escaping behavior.

(A) The trajectory of the worm for each scaling factor η. (B) Characteristics of the trajectory. The top graph represents $E_{θ} = {⟨ | θ_{i}^{(t)} - θ_{c t r l, i}^{(t)} | ⟩}_{i, t}$ . The middle graph shows the total traveled distance of the worm. The bottom graph represents the total absolute angle change ( $S = \sum_{t = 0}^{T - Δ t} | {⟨ s_{i}^{(t + Δ t)} ⟩}_{i} - {⟨ s_{i}^{(t)} ⟩}_{i} |$ , where T is the total time of the experimental video).

Appendix 1—figure 4

Download asset Open asset

The magnitude of the frictional force $| F_{b, i}^{(t)} |$ during the escaping behavior depending on the scaling factor η of the friction coefficient.

The color of each point in the heatmap represents the value of $| F_{b, i}^{(t)} |$ .

Appendix 1—figure 5

Download asset Open asset

The impact of the scaling factor η of the friction coefficients on the delta-turn.

(A) The trajectory of the worm for each scaling factor η. (B) Characteristics of the trajectory. The top graph represents $E_{θ} = {⟨ | θ_{i}^{(t)} - θ_{c t r l, i}^{(t)} | ⟩}_{i, t}$ . The middle graph shows the total traveled distance of the worm. The bottom graph represents the total absolute angle change ( $S = \sum_{t = 0}^{T - Δ t} | {⟨ s_{i}^{(t + Δ t)} ⟩}_{i} - {⟨ s_{i}^{(t)} ⟩}_{i} |$ )

Appendix 1—figure 6

Download asset Open asset

The magnitude of the frictional force $| F_{b, i}^{(t)} |$ during the delta-turn depending on the scaling factor η of the friction coefficient.

The color of each point in the heatmap represents the value of $| F_{b, i}^{(t)} |$ .

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. The source code of the program used in this study can be obtained from the open database GitHub (https://github.com/taegonchung/elegansbot, copy archived at Chung, 2024) or Python software repository PyPI (https://pypi.org/project/ElegansBot/).

The following previously published data sets were used

(2016) Dryad Digital Repository
Data from: Resolving coiled shapes reveals new reorientation behaviors in C. elegans.

https://doi.org/10.5061/dryad.t0m6p

References

1. Berri S
2. Boyle JH
3. Tassieri M
4. Hope IA
5. Cohen N
(2009) Forward locomotion of the nematode C. elegans is achieved through modulation of a single gait
HFSP Journal 3:186–193.

https://doi.org/10.2976/1.3082260
- PubMed
- Google Scholar
1. Boyle JH
2. Berri S
3. Cohen N
(2012) Gait modulation in C. elegans: an integrated neuromechanical model
Frontiers in Computational Neuroscience 6:10.

https://doi.org/10.3389/fncom.2012.00010
- PubMed
- Google Scholar
(2016) Resolving coiled shapes reveals new reorientation behaviors in C. elegans
eLife 5:e17227.

https://doi.org/10.7554/eLife.17227
- PubMed
- Google Scholar
Book
1. Butcher JC
(2003) Numerical Methods for Ordinary Differential Equations
John Wiley & Sons.

https://doi.org/10.1002/0470868279
- Google Scholar
Software
1. Chung T
(2024) Elegansbot, version swh:1:rev:2788144ad8194c9cbba726f445babd88dd6c00c0
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:6dda3ea0a880f9ec5ea53f72bf4795c227ee6366;origin=https://github.com/TaegonChung/elegansbot;visit=swh:1:snp:6363fbefb3688a04b797f3f52b2d6f985d7dccfc;anchor=swh:1:rev:2788144ad8194c9cbba726f445babd88dd6c00c0
1. Cohen E
2. Yemini E
3. Schafer W
4. Feitelson DG
5. Treinin M
(2012) Locomotion analysis identifies roles of mechanosensory neurons in governing locomotion dynamics of C. elegans
The Journal of Experimental Biology 215:3639–3648.

https://doi.org/10.1242/jeb.075416
- PubMed
- Google Scholar
1. Cook SJ
2. Jarrell TA
3. Brittin CA
4. Wang Y
5. Bloniarz AE
6. Yakovlev MA
7. Nguyen KCQ
8. Tang LTH
9. Bayer EA
10. Duerr JS
11. Bülow HE
12. Hobert O
13. Hall DH
14. Emmons SW
(2019) Whole-animal connectomes of both Caenorhabditis elegans sexes
Nature 571:63–71.

https://doi.org/10.1038/s41586-019-1352-7
- PubMed
- Google Scholar
1. Ekeberg Ö
(1993) A combined neuronal and mechanical model of fish swimming
Biological Cybernetics 69:363–374.

https://doi.org/10.1007/BF00199436
- Google Scholar
Book
1. Hall DH
2. Altun ZF
(2008)
C Elegans Atlas

Cold Spring Harbor Laboratory Press.
- Google Scholar
1. Harris CR
2. Millman KJ
3. van der Walt SJ
4. Gommers R
5. Virtanen P
6. Cournapeau D
7. Wieser E
8. Taylor J
9. Berg S
10. Smith NJ
11. Kern R
12. Picus M
13. Hoyer S
14. van Kerkwijk MH
15. Brett M
16. Haldane A
17. Del Río JF
18. Wiebe M
19. Peterson P
20. Gérard-Marchant P
21. Sheppard K
22. Reddy T
23. Weckesser W
24. Abbasi H
25. Gohlke C
26. Oliphant TE
(2020) Array programming with numpy
Nature 585:357–362.

https://doi.org/10.1038/s41586-020-2649-2
- PubMed
- Google Scholar
(2021) WormPose: Image synthesis and convolutional networks for pose estimation in C. elegans
PLOS Computational Biology 17:e1008914.

https://doi.org/10.1371/journal.pcbi.1008914
- PubMed
- Google Scholar
1. Hill AV
(1938) The heat of shortening and the dynamic constants of muscle
Proceedings of the Royal Society of London. Series B - Biological Sciences 126:136–195.

https://doi.org/10.1098/rspb.1938.0050
- Google Scholar
1. Hu DL
2. Nirody J
3. Scott T
4. Shelley MJ
(2009) The mechanics of slithering locomotion
PNAS 106:10081–10085.

https://doi.org/10.1073/pnas.0812533106
- PubMed
- Google Scholar
1. Hunter JD
(2007) Matplotlib: a 2d graphics environment
Computing in Science & Engineering 9:90–95.

https://doi.org/10.1109/MCSE.2007.55
- Google Scholar
1. Izquierdo EJ
2. Beer RD
(2018) From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegans
Philosophical Transactions of the Royal Society B: Biological Sciences 373:20170374.

https://doi.org/10.1098/rstb.2017.0374
- Google Scholar
1. Javer A
2. Currie M
3. Lee CW
4. Hokanson J
5. Li K
6. Martineau CN
7. Yemini E
8. Grundy LJ
9. Li C
10. Ch’ng Q
11. Schafer WR
12. Nollen EAA
13. Kerr R
14. Brown AEX
(2018) An open-source platform for analyzing and sharing worm-behavior data
Nature Methods 15:645–646.

https://doi.org/10.1038/s41592-018-0112-1
- PubMed
- Google Scholar
1. Jung J
2. Nakajima M
3. Takeuchi M
4. Najdovski Z
5. Huang Q
6. Fukuda T
(2016) Microfluidic device to measure the speed of C. elegans using the resistance change of the flexible electrode
Micromachines 7:50.

https://doi.org/10.3390/mi7030050
- PubMed
- Google Scholar
Conference
(2015) Numba: A LLVM-based Python JIT compiler
Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC. pp. 1–6.

https://doi.org/10.1145/2833157.2833162
- Google Scholar
1. Niebur E
2. Erdös P
(1991) Theory of the locomotion of nematodes: dynamics of undulatory progression on a surface
Biophysical Journal 60:1132–1146.

https://doi.org/10.1016/S0006-3495(91)82149-X
- PubMed
- Google Scholar
(2012) Dopamine signaling is essential for precise rates of locomotion by C. elegans
PLOS ONE 7:e38649.

https://doi.org/10.1371/journal.pone.0038649
- PubMed
- Google Scholar
(2013) Shifts in the Distribution of Mass Densities Is a Signature of Caloric Restriction in Caenorhabditis elegans
PLOS ONE 8:e69651.

https://doi.org/10.1371/journal.pone.0069651
- Google Scholar
1. Sakamoto K
2. Soh Z
3. Suzuki M
4. Iino Y
5. Tsuji T
(2021) Forward and backward locomotion patterns in C. elegans generated by a connectome-based model simulation
Scientific Reports 11:13737.

https://doi.org/10.1038/s41598-021-92690-2
- PubMed
- Google Scholar
1. Savitzky A
2. Golay MJE
(1964) Smoothing and differentiation of data by simplified least squares procedures
Analytical Chemistry 36:1627–1639.

https://doi.org/10.1021/ac60214a047
- Google Scholar
(2012) Undulatory locomotion of Caenorhabditis elegans on wet surfaces
Biophysical Journal 102:2772–2781.

https://doi.org/10.1016/j.bpj.2012.05.012
- PubMed
- Google Scholar
(1995)
“Prosthetic condition” and satisfaction with dentures

Nederlands Tijdschrift Voor Tandheelkunde 102:12–14.
- PubMed
- Google Scholar
1. Vidal-Gadea A
2. Topper S
3. Young L
4. Crisp A
5. Kressin L
6. Elbel E
7. Maples T
8. Brauner M
9. Erbguth K
10. Axelrod A
11. Gottschalk A
12. Siegel D
13. Pierce-Shimomura JT
(2011) Caenorhabditis elegans selects distinct crawling and swimming gaits via dopamine and serotonin
PNAS 108:17504–17509.

https://doi.org/10.1073/pnas.1108673108
- PubMed
- Google Scholar
1. Virtanen P
2. Gommers R
3. Oliphant TE
4. Haberland M
5. Reddy T
6. Cournapeau D
7. Burovski E
8. Peterson P
9. Weckesser W
10. Bright J
11. van der Walt SJ
12. Brett M
13. Wilson J
14. Millman KJ
15. Mayorov N
16. Nelson ARJ
17. Jones E
18. Kern R
19. Larson E
20. Carey CJ
21. Polat İ
22. Feng Y
23. Moore EW
24. VanderPlas J
25. Laxalde D
26. Perktold J
27. Cimrman R
28. Henriksen I
29. Quintero EA
30. Harris CR
31. Archibald AM
32. Ribeiro AH
33. Pedregosa F
34. van Mulbregt P
(2020) SciPy 1.0: fundamental algorithms for scientific computing in python
Nature Methods 17:261–272.

https://doi.org/10.1038/s41592-019-0686-2
- PubMed
- Google Scholar
(1986) The structure of the nervous system of the nematode Caenorhabditis elegans
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 314:1–340.

https://doi.org/10.1098/rstb.1986.0056
- PubMed
- Google Scholar
1. Yeon J
2. Kim J
3. Kim DY
4. Kim H
5. Kim J
6. Du EJ
7. Kang K
8. Lim HH
9. Moon D
10. Kim K
(2018) A sensory-motor neuron type mediates proprioceptive coordination of steering in C. elegans via two TRPC channels
PLOS Biology 16:e2004929.

https://doi.org/10.1371/journal.pbio.2004929
- PubMed
- Google Scholar

Article and author information

Author details

Taegon Chung

Daegu Gyeongbuk Institute of Science and Technology, Daegu, Republic of Korea

Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0004-2335-3500
Iksoo Chang

Daegu Gyeongbuk Institute of Science and Technology, Daegu, Republic of Korea

Contribution
Formal analysis, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing

Competing interests
No competing interests declared
Sangyeol Kim

Daegu Gyeongbuk Institute of Science and Technology, Daegu, Republic of Korea

Contribution
Conceptualization, Formal analysis, Validation, Investigation, Methodology, Writing - original draft, Writing - review and editing

For correspondence
sykim@dgist.ac.kr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0001-3200-9726

Funding

The p-CoE program of DGIST (23-CoE-BT-01)

Iksoo Chang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work is supported by the p-CoE program of DGIST 23-CoE-BT-01. We also acknowledge Prof. Kyuhyung Kim for fruitful discussions on experiments of C.elegans.

Version history

Sent for peer review: September 11, 2023
Preprint posted: September 27, 2023
Reviewed Preprint version 1: December 19, 2023
Reviewed Preprint version 2: March 27, 2024
Version of Record published: April 29, 2024

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.92562. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Taegon Chung
Iksoo Chang
Sangyeol Kim

(2024)

Development of equation of motion deciphering locomotion including omega turns of Caenorhabditis elegans

eLife 12:RP92562.

https://doi.org/10.7554/eLife.92562.3

Categories and tags

Research organism

C. elegans

Share this article

Cite this article

Components of ElegansBot.

Simulated locomotion from a sine kymogram.

Simulated locomotion from a kymogram of a real worm locomotion video.

Frictional force on each rod.

Body shape transition from the shape in water to the shape in agar.

Method of compressing motion state information.

The i-actuator.

The effect of the scaling factor η of the friction coefficients on the escaping behavior.

The magnitude of the frictional force |Fb,i(t)| during the escaping behavior depending on the scaling factor η of the friction coefficient.

The impact of the scaling factor η of the friction coefficients on the delta-turn.

The magnitude of the frictional force |Fb,i(t)| during the delta-turn depending on the scaling factor η of the friction coefficient.

Author details

Taegon Chung

Contribution

Competing interests

Iksoo Chang

Contribution

Competing interests

Sangyeol Kim

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Research organism

Further reading

The magnitude of the frictional force $| F_{b, i}^{(t)} |$ during the escaping behavior depending on the scaling factor η of the friction coefficient.

The magnitude of the frictional force $| F_{b, i}^{(t)} |$ during the delta-turn depending on the scaling factor η of the friction coefficient.