Introduction

With only a few hundred neurons, Caenorhabditis Elegans(C. elegans) performs various behaviors such as locomotion, sleeping, reproduction, and hunting (Hall & 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 & 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(Boyle, 2009; Boyle et al., 2012; Hu et al., 2009; Niebur & 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.

Kinetic simulators of rod-shaped animals such as C. elegans (Boyle et al., 2012; Niebur & 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 & 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 kinetic 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 kinetic 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 & Beer, 2018).

Even though there were studies on kinetic 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 kinetic simulators. A simulator should have an operational structure that imitates physical quantities from an experimental video to reproduce the motion of C. elegans. 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 2-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 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 (Fig. 1A). Let M (=2 μg) 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 (Fig. 1A). The mass, length, and moment of inertia of each rod is m = M/n, 2r = L/n, and I = mr²/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 the “i-rod”. The point where i-rod and (i+1)-rod meets is “i-joint”.

Figure 1.

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) (Fig. 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 Supplementary Information).

Value of time-dependent variable such as d_c, v_c, s_i, and ω_i at given time, t will be expressed as *^(t). When initial values, , and are given, the Newtonian equation of motion for acceleration, a_c and angular acceleration, {α_i}_{i∈{1,…,n}} must be acquired and numerically integrated twice to find , and 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 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 Supplementary Information). Because of this preservation of linearity, the frictional forces of translational motion (Fig. 1D) and rotational motion (Fig. 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 (Boyle, 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 (“·”: dot product of vectors, : unit vector parallel to i-rod, : unit vector perpendicular to i-rod), and the total frictional torque that i-rod receives is (positive or negative values are for torque pointing away from or into the paper plane, respectively.) (The proof is in Supplementary Information).

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 & 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 (Fig. 1C) so that there are 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, 1937). 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 θ_ctrl,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 (θ_ctrl,i) is a variable inside the elastic part of the muscle torque (τ_κ,i), to which τ_κ,i drives θ_i close. Also, control angle (θ_ctrl,i), which can be expressed by a heatmap (Fig. 2A, 2C, 3A, 3B), 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 Supplementary Information).

Figure 2.

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 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 Supplementary Information, Fig. S2).

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 (Fig. 1F) which we name as “joint force”. The joint force that (i+1)-rod exerts on i-rod is symbolized as F_i ≡ 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 (Fig. 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 Supplementary Information). When F₀ = F_n = 0, then the total joint force that i-rod receives is F_joint,i = F_i − F_i−1 and the total torque caused by joint force is τ_joint,i = [r_i × (F_i + F_i+1)] where “×” between two vectors means cross product.

As all forces and torques are found, 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⁻⁵ sec. Because the only external force exerts on the worm is the frictional force, the equation of translational motion is . If friction coefficients (b, b) are significantly greater than , numerical integration using explicit Euler method becomes unstable (Butcher, 2004). So, we tackled this instability of numerical integration by developing semi-implicit Euler Method , which makes numerical integration stable when any frictional coefficients greater than or equal to 0 is given (Details in Supplementary Information).

The equation of rotational motion of i-rod is Iα_i = τ_total,i = τ_b,i + τ_cκ,i + τ_joint,i. When the friction-related value (b_∥r²), elasticity-related value (κΔt), or damping coefficient(c) is significantly larger than , numerical integration using explicit Euler method becomes unstable (Butcher, 2004). To solve this instability, we constructed a semi-implicit Euler method for rotational motion and an error-corrected equation for angular momentum (Details in Supplementary Information). By using these semi-implicit Euler methods, solutions for d of a worm at a given time can be available for the ground surface of agar whose b_⊥, b_∥ are significantly larger than , 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, 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 (v) and the period (T) are, respectively, 1.832 and 1.6 (sec) on the agar plate and 0.667 and 0.4 (sec) 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 (Fig. 2A) and swimming (Fig. 2C) was calculated by substituting amplitude (A), wavenumber (v), and period (T) into .

Crawling trajectory, which performs sinusoidal locomotion in the positive x-axis direction, was obtained by inputting a crawling kymogram as input to ElegansBot (Fig. 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 (Fig. 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 (Fig. S3A). 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/sec), 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 (Fig. 2D). The worm gained significant momentum in the forward direction of locomotion when the body bent in the c-shape (Fig. S3B). 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 (Fig. S3B), and the average velocity was 0.223 (mm/sec).

ElegansBot exhibits more complex behavior including the turn motion

Unlike previous C. elegans body kinetic simulation studies, our simulation can replicate the worm’s behavior using a kymogram (Fig. 3A, 3B) 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 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 paper (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.

Figure 3.

The trajectory (Fig. 3C, 3D, Video S1, S2) 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 (Fig. 3C, 3D). 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.

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 kinetic 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_water,⊥ = 5.2 × 10³ (μg/sec) and b_water,∥ = b_water,⊥/1.5, respectively, while in agar, these values are b_agar,⊥ = 1.28 × 10⁸ (μg/sec) and b_agar,∥ = b_agar,⊥/40 (Boyle et al., 2012). For environmental index σ ∈ [0,1], we have defined the vertical and horizontal friction coefficients in the environment between water (σ = 0) and agar (σ = 1) as and , respectively.

Under an environmental index σ, for various pairs of frequency-period (v, T) when the control angle is (with A = 0.6 (rad)), we have found the (v, T) that maximizes the worm’s average speed(optimal (v, T)) (Fig. 4A). The optimal (v, T) exhibits nearly linear distribution (Fig. 4B). We noticed a transition from swimming body shape to crawling body shape as σ varies (Fig. 4, S4). The optimal (v, T) for σ=0(water) is (0.65, 0.4 sec), matching the actual (v, T) value of swimming behavior (Vidal-Gadea et al., 2011). The optimal (v, T) for σ=1(agar) is (1.9, 0.8 sec), and the optimal v (1.9) matches the actual v value (1.832) for crawling behavior (Vidal-Gadea et al., 2011), with the optimal T (0.8 sec) being half the actual T value (1.6 sec).

Figure 4.

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 to 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 (Boyle, 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 & Beer, 2018; Niebur & 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

1. 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.

2. 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.

3. 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 seconds), the body angles for that period was predicted using linear interpolation.

4. 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 Rossum, 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. 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 (https://github.com/taegonchung/elegansbot), and the demo can be found at the laboratory website (https://biophysics.dgist.ac.kr/elegansbot/demo). This code calculated the ElegansBot simulation of 10 seconds of simulation time in about 10 seconds of run-time on an Intel E3-1230v5 CPU.

5. 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 Supplementary Information).

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.