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.sa0Useful: Findings that have focused importance and scope
- Landmark
- Fundamental
- Important
- Valuable
- Useful
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
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 (=2 µg, details in ‘Worm’s mass, actuator elasticity coefficient, and damping coefficient’ of Appendix) be the mass and (=1 mm) be the length of the worm. Midline was approximated as (=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 , , and , 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.’.
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 (), velocity vector (), the angle measured counterclockwise from the positive x-axis to the tangential direction of the rod () (Figure 1B), and angular velocity () of i-rod at a given time . However, the minimum information required to describe the motion of all rods includes the displacement vector () and velocity vector () of the worm’s center of mass, and for each rod (Details in ‘Minimum information required to describe the motion of each rod’ of Appendix).
Value of time-dependent variables such as , , , and at a given time, will be expressed as . When initial values, , , , and are given, the Newtonian equation of motion for acceleration, and angular acceleration, must be acquired and numerically integrated twice to find , , , and at a given time, . 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 and for a straightened worm, respectively. Each rod experiences of the frictional force the worm gets. Thus, the perpendicular and parallel friction coefficients of each rod are , , respectively. The ratio of perpendicular friction coefficient to parallel friction coefficient () 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 (‘·’: 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 ‘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 () 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() 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 that the elastic term is and the damping term is where is control angle, , and and are the elasticity and damping coefficients of an actuator, respectively.
Control angle () is a variable inside the elastic part of the muscle torque (), to which drives close. Also, the control angle (), 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. represents the damping effect of muscle cells and somatic cells near i-actuator. The elasticity coefficient () and damping coefficient () 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).
By assuming that i-rod receives torque () from i-actuator, and (i+1)-rod receives torque (), 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 . The total muscle force () that i-rod receives from both of its ends is as follows (Details in ‘Proof of muscle force’ of Appendix).
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 . By Newton’s third law of motion about action and reaction, the joint force that i-rod exerts on (i+1)-rod is (Figure 1F). Joint force () can be calculated from the previously introduced given values ( , , , , ) (Details in ‘Joint force calculation method’ of Appendix). When , then the total joint force that i-rod receives is and the total torque caused by joint force is 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 () used in this work is s unless otherwise noted. Because the only external force exerts on the worm is the frictional force, the equation of translational motion is . If friction coefficients ( , ) are significantly greater than , numerical integration using the explicit Euler method () becomes unstable (Butcher, 2003). 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 ‘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 . When the friction-related value (), elasticity-related value (), or damping coefficient() is significantly larger than , numerical integration using explicit Euler method () 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 , , , of a worm at a given time can be available for the ground surface of agar whose , 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, . 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 () 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 () 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 (), wavenumber (), and period () into into .
Crawling trajectory, which performs sinusoidal locomotion in the positive x-axis direction, was obtained by inputting a crawling kymogram as 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 () 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 and 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 and 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 as . 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 ( where i=4 to 15) was 3536 pN, compared to the average magnitude of the remaining parts ( = 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 ( where i=16 to 25) was 10,497 pN, compared to the average magnitude of the remaining parts ( = 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.
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 (μg/sec) and , respectively, while in agar, these values are (μg/sec) and (Boyle et al., 2012). For environmental index , we have defined the vertical and horizontal friction coefficients in the environment between water () and agar () as and , respectively.
Under an environmental index , for various pairs of frequency-period (, ) when the control angle is (with = 0.6 (rad)), we have found the (, ) that maximizes the worm’s average velocity(optimal (, )) (Figure 5A). The optimal (, ) 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 (, ) for =0(water) is (0.65, 0.4 s), matching the actual (, ) value of swimming behavior (Vidal-Gadea et al., 2011). The optimal (, ) 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).
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 and time as . When the escaping behavior kymogram input 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 decreased as decreased (Figure 5—figure supplement 2B). From = 1.0 to = 0.1, the total absolute angular change ( where is the total time of the experimental video.) increased as decreased. However, from = 0.7 to = 0, remained constant within the error of 0.33 rad, and the total traveled distance() decreased as decreased. The maximum total traveled distance was at = 0.8. Using the same analysis method with the kymogram input 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 also decreased as decreased (Figure 5—figure supplement 3B). From = 1 to = 0.6, increased as decreased. From = 0.6 to = 0, 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.e., ). We defined a sine function fitting for the body angle as (where and ). Let us denote the set of a quantity for all as . For a given time , by curve fitting the function to the set of body angles , we can obtain the parameters ( , , , ) (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 ( , , , ) as follows. For , we calculated ( , , , ) using the following equations:
Using these initial guess values, we curve-fitted to to obtain ( , , , ). For , we obtained the initial guess values for curve fitting as , , , . Then, using these initial guess values, we curve-fitted to to obtain ( , , , ). Since the phase is not continuous for all time , we defined a continuous value for all time as follows (Figure 4—figure supplement 1):
To obtain the derivatives of the noise-reduced smoothed values and for the raw data and , 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 and . We then calculated the temporal integrals and . Let us denote the average of a quantity for all time t as . We calculated and (Figure 4—figure supplement 1). Finally, we defined the worm’s behavior classification as turn when , as forward locomotion when and , and as backward locomotion when and (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 . Assuming that the worm is a cylindrical body with a bottom surface radius of 25µm, the volume of the worm is , and the density of the worm is (Reina et al., 2013), so the weight of the worm is .
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 , and the maximum value of the muscle damping coefficient is . When 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 is the torque elasticity coefficient , so . In ElegansBot, it was assumed that this value is constant regardless of , , . In the same way, the torque muscle damping coefficient is .
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 (), velocity (), angle () measured counterclockwise from the positive x-axis to the direction of vector , and angular velocity() of every i-rod. However, knowing only the position and velocity of the worm ( , ) and the angle and angular velocity of every i-rod ( , ) 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 , and the vector parallel to i-rod is . If the position vector where i-rod and (i+1)-rod meet is , and the free end of 1-rod and n-rod are and , then (Appendix 1—figure 1A).
3. Method of calculating the relative and absolute position of the rod
If the relative coordinate with as the origin is designated as , then the following equations satisfy.
Thus, can be calculated from and (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, . In the same way as the location information compression, the relative velocity vector with as the origin has the following relationships ().
Thus, can be calculated from , , and (Appendix 1—figure 1C).
Preservation of linearity in friction
The velocity of an arbitrary point particle which is included by i-rod can be decomposed into two velocity components . In this case, if the object is subject to an anisotropic Stokes friction, there is linearity between the friction , obtained from each velocity component , and the friction obtained from velocity .
Frictional torque by rotational motion
Let us denote variable as the distance from the center of i-rod measured along the direction of vector . It is to be noted that is within the range . For an infinitesimal where , the moment arm vector for the infinitesimal interval (hereafter referred to as the infinitesimal interval ) from the center of i-rod is . The coefficient of friction for the infinitesimal interval is:
The velocity component due to the rotational motion of the infinitesimal interval is . Therefore, the frictional force received by the infinitesimal interval due to rotational motion is:
The torque received by the infinitesimal interval due to rotational motion is:
The total frictional torque received by i-rod due to rotational motion is
Proof of muscle force
The muscle force, which makes the torque 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 ( , ) 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 (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).
The force that i-rod receives from i-actuator is and the magnitude of the force is and the direction of the force is . These equations come down to the following equations.
Therefore, the total force that i-rod receives from i-actuator and (i-1)-actuator is
Joint force calculation method
Let us calculate the joint force from the given values ( , , , , ). In this part, the superscript 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.
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.
Multiplying both sides by gives:
If is the total force applied to i-rod other than and , which is , then:
If is the total torque applied to i-rod excluding , which is , then:
Dividing both sides by gives:
If , then:
Following the same method,
To organize the terms of the equations into known values and unknown values,
Therefore, if we know all and for each i, we can find . If , and if we expand ,
If we set and , then:
If we set and express the above equation in a multidimensional tensor form,
Let us set:
As and can be calculated from , from , and from , we can find from .
Therefore, if we know , , , , we can find .
We can find , thus (), 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 , then:
And the matrix equation equivalent to this tensor equation is:
Let us set:
Since , are symmetric matrices, is a heptadiagonal symmetric matrix. Since 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 .
As a result, we can find the joint force from the known values ( , , , , ).
Proof of numerical integration for the translational motion of a worm using semi-implicit Euler method
When the friction coefficients , are sufficiently large compared to , numerical integration via the explicit Euler method () becomes unstable (Butcher, 2003). Therefore, for all friction coefficients greater than or equal to 0, the semi-implicit Euler method () 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.
The integration formula for using the implicit Euler method is as follows. (where )
Because and are unknown at time , the numerical calculation of the above formula is impossible.
If the above formula is assumed to be true, the following approximation can be used.
The approximate value of by the above approximation is as follows.
The approximate value of by the approximate value of is as follows.
If both sides are divided by ,
, satisfy the following because they are the internal forces of the worm ( is a zero vector).
Therefore, the approximation of the force received by the worm is as follows.
If both sides are divided by ,
This approximation ensures computational stability regardless of the size of , . That is, this approximation solves the problem of the decrease in computational stability of numerical integration through the explicit Euler method when , are sufficiently large compared to .
Numerical integration of the rotational motion of i-rod using semi-implicit Euler method
First, the numerical integration formulas for and using the implicit Euler method are as follows.
If we set , the equation describing the rotation of i-rod is as follows.
If the above formula is expanded for time ,
The above formula is impossible to integrate because , , are unknown at time t.
If the above formula is assumed to be true, the following approximation can be used.
The above formula can be expressed as a matrix formula as follows.
In the above vector matrix formula, let us represent the vectors and matrix by the following symbols.
Then, the matrix formula is expressed as follows.
Now, the following approximation can be obtained where is a unit matrix of size .
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 , the calculation error accumulated for the rotational inertia of the whole worm significantly influenced the calculation result and . To prevent this, the error is corrected as follows. If , , the moment of inertia of the entire worm at time t is as follows by the parallel axis theorem.
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)
If i-rod is approximated as a point particle, the rotational inertia of the whole worm at time is as follows.
The predicted value of at , , is calculated by the semi-implicit Euler method (See ‘Numerical integration of the rotational motion’ in Appendix). ()
The predicted value of at is as follows.
Predicted values , , for , , are calculated from , , , (See ‘Minimum information required to describe the motion of each rod’ in Appendix).
Using , , , the moment of inertia and rotational inertia of the entire worm at time are calculated as follows. (where , , )
is calculated as follows.
is calculated as follows.
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 , as if .
Proper selection of friction coefficients
When using the vertical and horizontal friction coefficients , 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 , by multiplying scaling factor to the agar friction coefficients , of the previous work (Boyle et al., 2012). Let us denote the set of a quantity for all pairs of index i and time as . When the kymogram input 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 () was set to 10-6 sec for higher accuracy of simulation). Let us denote the average of a quantity for all pairs as . When was 1 ∼ 10-9 , the smaller , the smaller was. The reduction in 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, did not fall below 0.125 (rad). This is because there was a time delay between the input and the response . As decreased from 1 to 10-2 , the total traveled distance of the worm () and the total absolute angle change () 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, , and the pattern of (Appendix 1—figure 4) were similar across trajectories. If was smaller than 10-6 , the worm’s total traveled distance decreased and 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 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 was (Appendix 1—figure 5B). The reduction in 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, did not fall below 0.15 (rad). As decreased from 1 to 10-2 , 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, , and the pattern of (Appendix 1—figure 6) were similar across trajectories. If was smaller than 10-6 , the worm’s total traveled distance decreased and 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(), selecting an appropriate value (between 10-2 and 10-6) was crucial for replicating the trajectory of sequenced locomotive behavior. Therefore, we chose , () as the friction coefficients that sufficiently reduce among those closest to the agar friction coefficients of the previous work (Boyle et al., 2012) for both escaping behavior and delta-turn.
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/).
-
Dryad Digital RepositoryData from: Resolving coiled shapes reveals new reorientation behaviors in C. elegans.https://doi.org/10.5061/dryad.t0m6p
References
-
Gait modulation in C. elegans: an integrated neuromechanical modelFrontiers in Computational Neuroscience 6:10.https://doi.org/10.3389/fncom.2012.00010
-
BookNumerical Methods for Ordinary Differential EquationsJohn Wiley & Sons.https://doi.org/10.1002/0470868279
-
Locomotion analysis identifies roles of mechanosensory neurons in governing locomotion dynamics of C. elegansThe Journal of Experimental Biology 215:3639–3648.https://doi.org/10.1242/jeb.075416
-
A combined neuronal and mechanical model of fish swimmingBiological Cybernetics 69:363–374.https://doi.org/10.1007/BF00199436
-
WormPose: Image synthesis and convolutional networks for pose estimation in C. elegansPLOS Computational Biology 17:e1008914.https://doi.org/10.1371/journal.pcbi.1008914
-
The heat of shortening and the dynamic constants of muscleProceedings of the Royal Society of London. Series B - Biological Sciences 126:136–195.https://doi.org/10.1098/rspb.1938.0050
-
The mechanics of slithering locomotionPNAS 106:10081–10085.https://doi.org/10.1073/pnas.0812533106
-
Matplotlib: a 2d graphics environmentComputing in Science & Engineering 9:90–95.https://doi.org/10.1109/MCSE.2007.55
-
From head to tail: a neuromechanical model of forward locomotion in Caenorhabditis elegansPhilosophical Transactions of the Royal Society B: Biological Sciences 373:20170374.https://doi.org/10.1098/rstb.2017.0374
-
ConferenceNumba: A LLVM-based Python JIT compilerProceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC. pp. 1–6.https://doi.org/10.1145/2833157.2833162
-
Smoothing and differentiation of data by simplified least squares proceduresAnalytical Chemistry 36:1627–1639.https://doi.org/10.1021/ac60214a047
-
Undulatory locomotion of Caenorhabditis elegans on wet surfacesBiophysical Journal 102:2772–2781.https://doi.org/10.1016/j.bpj.2012.05.012
-
“Prosthetic condition” and satisfaction with denturesNederlands Tijdschrift Voor Tandheelkunde 102:12–14.
-
The structure of the nervous system of the nematode Caenorhabditis elegansPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 314:1–340.https://doi.org/10.1098/rstb.1986.0056
Article and author information
Author details
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:
- Preprint posted:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
- Version of Record published:
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
© 2023, Chung et al.
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.
Metrics
-
- 711
- views
-
- 74
- downloads
-
- 0
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading
-
- Computational and Systems Biology
- Neuroscience
The basolateral amygdala (BLA) is a key site where fear learning takes place through synaptic plasticity. Rodent research shows prominent low theta (~3–6 Hz), high theta (~6–12 Hz), and gamma (>30 Hz) rhythms in the BLA local field potential recordings. However, it is not understood what role these rhythms play in supporting the plasticity. Here, we create a biophysically detailed model of the BLA circuit to show that several classes of interneurons (PV, SOM, and VIP) in the BLA can be critically involved in producing the rhythms; these rhythms promote the formation of a dedicated fear circuit shaped through spike-timing-dependent plasticity. Each class of interneurons is necessary for the plasticity. We find that the low theta rhythm is a biomarker of successful fear conditioning. The model makes use of interneurons commonly found in the cortex and, hence, may apply to a wide variety of associative learning situations.
-
- Cancer Biology
- Computational and Systems Biology
Effects from aging in single cells are heterogenous, whereas at the organ- and tissue-levels aging phenotypes tend to appear as stereotypical changes. The mammary epithelium is a bilayer of two major phenotypically and functionally distinct cell lineages: luminal epithelial and myoepithelial cells. Mammary luminal epithelia exhibit substantial stereotypical changes with age that merit attention because these cells are the putative cells-of-origin for breast cancers. We hypothesize that effects from aging that impinge upon maintenance of lineage fidelity increase susceptibility to cancer initiation. We generated and analyzed transcriptomes from primary luminal epithelial and myoepithelial cells from younger <30 (y)ears old and older >55y women. In addition to age-dependent directional changes in gene expression, we observed increased transcriptional variance with age that contributed to genome-wide loss of lineage fidelity. Age-dependent variant responses were common to both lineages, whereas directional changes were almost exclusively detected in luminal epithelia and involved altered regulation of chromatin and genome organizers such as SATB1. Epithelial expression of gap junction protein GJB6 increased with age, and modulation of GJB6 expression in heterochronous co-cultures revealed that it provided a communication conduit from myoepithelial cells that drove directional change in luminal cells. Age-dependent luminal transcriptomes comprised a prominent signal that could be detected in bulk tissue during aging and transition into cancers. A machine learning classifier based on luminal-specific aging distinguished normal from cancer tissue and was highly predictive of breast cancer subtype. We speculate that luminal epithelia are the ultimate site of integration of the variant responses to aging in their surrounding tissue, and that their emergent phenotype both endows cells with the ability to become cancer-cells-of-origin and represents a biosensor that presages cancer susceptibility.