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

  1. Taegon Chung
  2. Iksoo Chang
  3. Sangyeol Kim  Is a corresponding author
  1. Daegu Gyeongbuk Institute of Science and Technology, Republic of Korea

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

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 , 2r=L/n , and I=mr2/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.’.

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 (di), velocity vector (vi), the angle measured counterclockwise from the positive x-axis to the tangential direction of the rod (si) (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 (dc) and velocity vector (vc) of the worm’s center of mass, si 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 dc , vc , si , and ωi at a given time, t will be expressed as (t). When initial values, dc(0) , vc(0) , si(0) , and ωi(0) are given, the Newtonian equation of motion for acceleration, ac and angular acceleration, {αi}i1,,n must be acquired and numerically integrated twice to find dc(t) , vc(t) , si(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 Fb,i=bn(vir^i)r^ibn(viN^i)N^i (‘·’: dot product of vectors, r^i[cos(si)sin(si)]T : unit vector parallel to i-rod, N^i[sin(si)cos(si)]T : unit vector perpendicular to i-rod), and the total frictional torque that i-rod receives is τb,i=13bnr2ω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 θctrl,i is control angle, θi=si+1-si , 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, the control angle (θctrl,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
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
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 (Fb,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 (Fcκ,i) that i-rod receives from both of its ends is as follows (Details in ‘Proof of muscle force’ of Appendix).

Fcκ,i=j=10(1)j1τi+jsinθi+j2rcos2θi+j2[cos(si+j+si+j+12)sin(si+j+si+j+12)]

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 FiF(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 Fi(i+1)=F(i+1)i=Fi (Figure 1F). Joint force (Fi) can be calculated from the previously introduced given values (si , Fcκ,i , Fb,i , τcκ,i , τb,i) (Details in ‘Joint force calculation method’ of Appendix). When F0=Fn=0 , then the total joint force that i-rod receives is Fjoint,i=FiFi1 and the total torque caused by joint force is τjoint,i=[ri×(Fi+Fi+1)] where ‘×’ between two vectors means cross-product.

As all forces and torques are found, dc(t) , vc(t) , si(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 ac=iFb,iM . If friction coefficients (b , b) are significantly greater than MΔt , numerical integration using the explicit Euler method (vc(t+Δt)=vc(t)+ac(t)Δt=vc(t)+iFb,iMΔt) becomes unstable (Butcher, 2003). So, we tackled this instability of numerical integration by developing semi-implicit Euler Method (vc(t+Δt)=vc(t)+ac(t+Δt)Δtvc(t)+11+bΔtMiFb,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=τtotal,i=τb,i+τcκ,i+τjoint,i . When the friction-related value (br2), elasticity-related value (κΔt), or damping coefficient(c) is significantly larger than IΔt , numerical integration using explicit Euler method (ωi(t+Δt)=ωi(t)+αi(t)Δt=ωi(t)+τtotal,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 dc(t) , vc(t) , si(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 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{1,,n1}) 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 θctrl,i(t)=Acos(2π(ν(i1)/(n2)t/T)) .

Crawling trajectory, which performs sinusoidal locomotion in the positive x-axis direction, was obtained by inputting a crawling kymogram as θctrl,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 (θctrl,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 (|Fb,i(t)|i where i=4 to 15) was 3536 pN, compared to the average magnitude of the remaining parts (|Fb,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 (|Fb,i(t)|i where i=16 to 25) was 10,497 pN, compared to the average magnitude of the remaining parts (|Fb,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
Frictional force on each rod.

(A) Escaping behavior. The top panel represents the frictional force Fb,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 |Fb,i(t)| . The third panel from the top represents the average |Fb,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 bwater,=5.2×103 (μg/sec) and bwater,=bwater,/1.5 , respectively, while in agar, these values are bagar,=1.28×108(μg/sec) and bagar,=bagar,/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 bσ,=bwater,1σbagar,σ and bσ,=bwater,1σbagar,σ , respectively.

Under an environmental index σ, for various pairs of frequency-period (ν, T) when the control angle is θctrl,i(t)=Acos(2π(νi123tT)) (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
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 {θctrl,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)θctrl,i(t)|i,t decreased as σ decreased (Figure 5—figure supplement 2B). From σ = 1.0 to σ = 0.1, the total absolute angular change (S=t=0TΔt|si(t+Δt)isi(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(t|vc(t)Δt|) decreased as σ decreased. The maximum total traveled distance was at σ = 0.8. Using the same analysis method with the kymogram input {θctrl,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(i1)/(n2) (i.e., 0ξi1). We defined a sine function fitting for the body angle θi(t) as θ^i(t)=A(t)sin((2π/λ(t))(ξiξ0(t)))+θ0(t) (where A(t)0 and 0ξ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 θ^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 (A^(t) , λ^(t) , ξ^0(t) , θ^0(t)) as follows. For t=0, we calculated (A^(0) , λ^(0) , ξ^0(0) , θ^0(0)) using the following equations:

A^(0)=12(maxi(θi(0))mini(θi(0)))λ^(0)=2n2(argmaxi(θi(0))argmini(θi(0)))ξ^0(0)=1n2(argmaxi(θi(0)))λ^(0)4θ^0(0)=1n1(iθi(0))

Using these initial guess values, we curve-fitted θ^i(0) to {θi(0)}i to obtain (A(0) , λ(0) , ξ0(0) , θ0(0)). For tΔt, we obtained the initial guess values for curve fitting as A^(t)=A(tΔt) , λ^(t)=λ(tΔt) , ξ^0(t)=ξ0(tΔt) , θ^0(t)=θ0(tΔt) . Then, using these initial guess values, we curve-fitted θ^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 ξ~0(t) for all time t as follows (Figure 4—figure supplement 1):

ξ~0(t)={ξ0(t)if t=0ξ0(t)ξ0(tΔt)+ξ~0(tΔt)if t>0 and λ(t)2ξ0(t)ξ0(tΔt)λ(t)2ξ0(t)ξ0(tΔt)+ξ~0(tΔt)+λ(t)if t>0 and ξ0(t)ξ0(tΔt)<λ(t)2ξ0(t)ξ0(tΔt)+ξ~0(tΔt)λ(t)if t>0 and λ(t)2<ξ0(t)ξ0(tΔt)

To obtain the derivatives of the noise-reduced smoothed values ξ¯0(t) and θ¯0(t) for the raw data ξ~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 dξ¯0(t)dt and dθ¯0(t)dt . We then calculated the temporal integrals ξ0(ζ)t=0ζdξ¯0(t)dtΔt and θ0(ζ)t=0ζdθ¯0(t)dtΔt. Let us denote the average of a quantity for all time t as t. We calculated ξ¯0(t)=ξ0(t)ξ0(t)t+ξ~0(t)t and θ¯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 θ¯0(t)<0.07, as forward locomotion when θ¯0(t)0.07 and dξ¯0(t)dt>0, and as backward locomotion when θ¯0(t)0.07 and dξ¯0(t)dt0 (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×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 1mm0.025mm2π0.002mm3 , and the density of the worm is 1000μg/mm3 (Reina et al., 2013), so the weight of the worm is M=1000μg/mm3×0.002mm3=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 kmax=2.8108μg/sec2 , and the maximum value of the muscle damping coefficient is kmax/5.61sec . 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 κ=dτκ,idθi=ddθi(γ(k(2γtan(θi/2))))kγ2=1.75105μgmm2/(sec2rad) . 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.61.75105μgmm2/secrad .

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 (di), velocity (vi), angle (si) measured counterclockwise from the positive x-axis to the direction of vector ri , and angular velocity(ωi) of every i-rod. However, knowing only the position and velocity of the worm (dc , vc) and the angle and angular velocity of every i-rod (si , ω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 dc=[xcyc]T=(i=1nmdi)/M=(i=1ndi)/n, and the vector parallel to i-rod is rir[cos(si)sin(si)]T=rr^i . If the position vector where i-rod and (i+1)-rod meet is dtipi , and the free end of 1-rod and n-rod are dtip0 and dtipn , then dtipi=di+ri=di+1-ri+1 (Appendix 1—figure 1A).

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

If the relative coordinate with dtip0 as the origin is designated as d`, then the following equations satisfy.

dtip0=[00]dtipi=j=1i2rjdi=dtipi+dtipi12dc=1ni=1ndidi=didc+dc

Thus, di can be calculated from dc and si (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, vc=(i=1nmvi)/M=(i=1nvi)/n. In the same way as the location information compression, the relative velocity vector v` with vtip0 as the origin has the following relationships (z^r^i×N^i).

vtip0=[00]vtipi=j=1iωjz^×(2rj)vi=vtipi+vtipi12vc=1ni=1nvivi=vivc+vc

Thus, vi can be calculated from vc , si , 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 Fb,α , Fb,β obtained from each velocity component vα , vβ and the friction Fb obtained from velocity v.

Fb,α=bn1(vαr^i)r^ibn1(vαN^i)N^iFb,β=bn1(vβr^i)r^ibn1(vβN^i)N^iFb=bn1(vr^i)r^ibn1(vN^i)N^i=bn1((vα+vβ)r^i)r^ibn1((vα+vβ)N^i)N^i=bn1(vαr^i)r^ibn1(vαN^i)N^ibn1(vβr^i)r^ibn1(vβN^i)N^i=Fb,α+Fb,β

Frictional torque by rotational motion

Let us denote variable ρ as the distance from the center of i-rod measured along the direction of vector 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 ρr^i . The coefficient of friction for the infinitesimal interval ρ is:

bndρ2r

The velocity component due to the rotational motion of the infinitesimal interval ρ is ρωiN^i . Therefore, the frictional force received by the infinitesimal interval ρ due to rotational motion is:

bndρ2rρωiN^i=12bnrωiρdρN^i

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

dτ=(ρr^i)×(12bnrωiρdρN^i)=12bnrωiρ2dρr^i×N^i=12bnrωiρ2dρz^

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

τb,i=rrdτ=rr12bnrωiρ2dρz^=12bnrωi[13ρ3]rrz^=12bnrωi23r3z^=13bnr2ωiz^

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 (Fsup1,i , Fsup2,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).

Fsup1,i=Fp1,icosθi2τi=(ri×Fsup1,i)z^=rFsup1,icosθi2Fsup1,i=τircosθi2Fp1,i=Fsup1,icosθi2=τircos2θi2

The force that i-rod receives from i-actuator is Fp1,i+Fsup1,i and the magnitude of the force is Fp1,i+Fsup1,i=Fp1,isinθi/2=τisinθi/2/rcos2θi/2 and the direction of the force is [cos((si+si+1)/2)sin((si+si+1)/2)]T . These equations come down to the following equations.

Fp1,i+Fsup1,i=τisinθi2rcos2θi2[cos(si+si+12)sin(si+si+12)]
Fp1,i+Fsup1,i+Fp2,i+Fsup2,i=0
Fp2,i+Fsup2,i=τisinθi2rcos2θi2[cos(si+si+12)sin(si+si+12)]

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

Fcκ,i=j=10(1)j1τi+jsinθi+j2rcos2θi+j2[cos(si+j+si+j+12)sin(si+j+si+j+12)]

Joint force calculation method

Let us calculate the joint force Fi from the given values (si , Fcκ,i , Fb,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.

di+ri=di+1-ri+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.

d2dt2(di+ri)=d2dt2(di+rr^i)=ddt(vi+rωiN^i)=ai+rdωidtN^i+rωidN^idt=ai+rαiN^irωi2r^i=ai+αi×riωi2ri
d2dt2(di+1ri+1)=d2dt2(di+1rr^i+1)=ddt(vi+1rωi+1N^i+1)=ai+1rdωi+1dtN^i+1rωi+1dN^i+1dt=ai+1rαi+1N^i+1+rωi+12r^i+1=ai+1αi+1×ri+1+ωi+12ri+1
d2dt2(di+ri)=d2dt2(di+1ri+1)=ai+αi×riωi2ri=ai+1αi+1×ri+1+ωi+12ri+1

Multiplying both sides by m gives:

m(ai+αi×ri)mωi2ri=m(ai+1αi+1×ri+1)+mωi+12ri+1

If Fres,i is the total force applied to i-rod other than Fi and -Fi-1 , which is Fres,i=Fcκ,i+Fb,i , then:

mai=Fi-Fi-1+Fres,i

If τres,i is the total torque applied to i-rod excluding ri×Fi-1+Fi , which is τres,i=τcκ,i+τb,i , then:

Iαi=13mr2αi=ri×(Fi1+Fi)+τres,i

Dividing both sides by 13r2 gives:

mαi=3r2ri×(Fi1+Fi)+3r2τres,i

If hres,i3r2τres,i×ri , then:

m(ai+αi×ri)=(FiFi1+Fres,i)+[3r2ri×(Fi1+Fi)+3r2τres,i]×ri=(FiFi1)+3[r^i×(Fi1+Fi)]×r^i+Fres,i+3r2τres,i×ri=(FiFi1)+3[(Fi1+Fi)N^i]N^i+Fres,i+hres,i

Following the same method,

m(ai+1αi+1×ri+1)=(Fi+1Fi)3[(Fi+Fi+1)N^i+1]N^i+1+Fres,i+1hres,i+1

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

kiFiFi1Pi3N^iN^iThi3[(Fi1+Fi)N^i]N^i=3N^iN^iT(Fi1+Fi)=Pi(Fi1+Fi)m(ai+αi×ri)=ki+hi+Fres,i+hres,im(ai+1αi+1×ri+1)=ki+1hi+1+Fres,i+1hres,i+1ki+hi+Fres,i+hres,imωi2ri=ki+1hi+1+Fres,i+1hres,i+1+mωi+12ri+1qiki+hiki+1+hi+1=Fres,ihres,i+Fres,i+1hres,i+1+mωi2ri+mωi+12ri+1

Therefore, if we know all Fres,i and τres,i for each i, we can find qi . If I[1001] , and if we expand qi ,

qiki+hiki+1+hi+1=ki+1+ki+hi+hi+1=Fi+1+Fi+FiFi1+Pi(Fi1+Fi)+Pi+1(Fi+Fi+1)=Fi1(PiI)+Fi(Pi+Pi+1+2I)+Fi+1(Pi+1I)

If we set AiPi-I and BiPi+Pi+1+2I, then:

qi=AiFi-1+BiFi+Ai+1Fi+1

If we set 0¯=[0000] and express the above equation in a multidimensional tensor form,

[B1A20¯0¯0¯0¯A2B2A30¯0¯0¯0¯A3B3A40¯0¯0¯0¯A4B4A50¯0¯0¯0¯A5B5An10¯0¯0¯0¯An1Bn1][F1F2F3F4Fn2Fn1]=[q1q2q3q4qn2qn1]

Let us set:

D[B1A20¯0¯0¯0¯A2B2A30¯0¯0¯0¯A3B3A40¯0¯0¯0¯A4B4A50¯0¯0¯0¯A5B5An10¯0¯0¯0¯An1Bn1]
F[F1F2F3F4Fn2Fn1]
Q[q1q2q3q4qn2qn1]

As Ai and Bi can be calculated from Pi , Pi from N^i , and N^i from si , we can find D from si .

qi=Fres,i+Fres,i+1hres,ihres,i+1+mωi2ri+mωi+12ri+1=Fres,i+Fres,i+13r2(τres,i×ri+τres,i+1×ri+1)+mωi2ri+mωi+12ri+1

Therefore, if we know Fcκ,i , Fb,i , τcκ,i , τb,i , we can find Q.

DF=Q

We can find F, thus Fi (i{1,,n1}), 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:

[[(B1)1,1(B1)1,2(B1)2,1(B1)2,2][(A2)1,1(A2)1,2(A2)2,1(A2)2,2][0000][(A2)1,1(A2)1,2(A2)2,1(A2)2,2][(B2)1,1(B2)1,2(B2)2,1(B2)2,2][(A3)1,1(A3)1,2(A3)2,1(A3)2,2][0000][(A3)1,1(A3)1,2(A3)2,1(A3)2,2][(B3)1,1(B3)1,2(B3)2,1(B3)2,2]][[(F1)1,1(F1)2,1][(F2)1,1(F2)2,1][(F3)1,1(F3)2,1]]=[[(q1)1,1(q1)2,1][(q2)1,1(q2)2,1][(q3)1,1(q3)2,1]]

And the matrix equation equivalent to this tensor equation is:

[(B1)1,1(B1)1,2(A2)1,1(A2)1,200(B1)2,1(B1)2,2(A2)2,1(A2)2,200(A2)1,1(A2)1,2(B2)1,1(B2)1,2(A3)1,1(A3)1,2(A2)2,1(A2)2,2(B2)2,1(B2)2,2(A3)2,1(A3)2,200(A3)1,1(A3)1,2(B3)1,1(B3)1,200(A3)2,1(A3)2,2(B3)2,1(B3)2,2][(F1)1,1(F1)2,1(F2)1,1(F2)2,1(F3)1,1(F3)2,1]=[(q1)1,1(q1)2,1(q2)1,1(q2)2,1(q3)1,1(q3)2,1]

Let us set:

D[(B1)1,1(B1)1,2(A2)1,1(A2)1,200(B1)2,1(B1)2,2(A2)2,1(A2)2,200(A2)1,1(A2)1,2(B2)1,1(B2)1,2(A3)1,1(A3)1,2(A2)2,1(A2)2,2(B2)2,1(B2)2,2(A3)2,1(A3)2,200(A3)1,1(A3)1,2(B3)1,1(B3)1,200(A3)2,1(A3)2,2(B3)2,1(B3)2,2]

Since Ai , Bi 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 Fi .

As a result, we can find the joint force Fi from the known values (si , Fcκ,i , Fb,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 (vc(t+Δt)=vc(t)+ac(t)Δt=vc(t)+iFb,i(t)MΔt) becomes unstable (Butcher, 2003). Therefore, for all friction coefficients greater than or equal to 0, the semi-implicit Euler method (vc(t+Δt)=vc(t)+ac(t+Δt)Δtvc(t)+11+bΔtMiFb,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.

Fi=mai=Fb,i+Fcκ,i+Fjoint,i=bnv,ibnv,i+Fcκ,i+Fjoint,i

The integration formula for ai using the implicit Euler method is as follows. (where bb)

vi(t+Δt)=vi(t)+ai(t+Δt)Δt=vi(t)Δtm(bnv,i(t+Δt)+bnv,i(t+Δt))+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))=vi(t)ΔtM(bv,i(t+Δt)+bv,i(t+Δt))+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))=vi(t)ΔtM(b(v,i(t+Δt)+v,i(t+Δt))+(bb)v,i(t+Δt))+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))=vi(t)ΔtM(bvi(t+Δt)+(bb)v,i(t+Δt))+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))(1+bΔtM)vi(t+Δt)=vi(t)ΔtM((bb)v,i(t+Δt))+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))=vi(t)+(bb)ΔtMv,i(t+Δt)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))

Because v,i(t+Δt) and v,i(t+Δt) are unknown at time t, the numerical calculation of the above formula is impossible.

|(bb)ΔtM(v,i(t+Δt)v,i(t))||vi(t)+(bb)ΔtMv,i(t+Δt)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))|0

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

(1+bΔtM)vi(t+Δt)=vi(t)+(bb)ΔtMv,i(t+Δt)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))vi(t)+(bb)ΔtMv,i(t)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))

The approximate value of vi(t+Δt) by the above approximation is as follows.

vi(t+Δt)[vi(t)+(bb)ΔtMv,i(t)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))]/(1+bΔtM)

The approximate value of ai(t+Δt)Δt by the approximate value of vi(t+Δt) is as follows.

ai(t+Δt)Δt=vi(t+Δt)vi(t)[vi(t)+(bb)ΔtMv,i(t)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))]/(1+bΔtM)vi(t)=[bΔtMvi(t)+(bb)ΔtMv,i(t)+Δtm(Fcκ,i(t+Δt)+Fjoint,i(t+Δt))]/(1+bΔtM)=Δtm[bnvi(t)+(bb)nv,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)]/(1+bΔtM)=Δtm[bn(vi(t)v,i(t))bnv,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)]/(1+bΔtM)=Δtm[bnv,i(t)bnv,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)]/(1+bΔtM)=Δtm[Fb,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)]/(1+bΔtM)

If both sides are divided by Δt,

ai(t+Δt)1(1+bΔtM)Fb,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)m

Fcκ,i(t+Δt) , Fjoint,i(t+Δt) satisfy the following because they are the internal forces of the worm (0 is a zero vector).

iFcκ,i(t+Δt)=0iFjoint,i(t+Δt)=0

Therefore, the approximation of the force Mac(t+Δt) received by the worm is as follows.

Mac(t+Δt)=imai(t+Δt)im1(1+bΔtM)Fb,i(t)+Fcκ,i(t+Δt)+Fjoint,i(t+Δt)m=1(1+bΔtM)(iFb,i(t)+iFcκ,i(t+Δt)+iFjoint,i(t+Δt))=1(1+bΔtM)iFb,i(t)

If both sides are divided by M,

ac(t+Δt)1(1+bΔtM)iFb,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.

si(t+Δt)=si(t)+ωi(t+Δt)Δt
ωi(t+Δt)=ωi(t)+αi(t+Δt)Δt

If we set β13bnr2 , the equation describing the rotation of i-rod is as follows.

Iαi=τb,i+τcκ,i+τjoint,i=βωi+c(ωi+12ωi+ωi1)+κ(θiθi1θctrl,i+θctrl,i1)+τjoint,i=βωi+c(ωi+12ωi+ωi1)+κ(si+12si+si1θctrl,i+θctrl,i1)+τjoint,i=cωi1(β+2c)ωi+cωi+1+κ(si12si+si+1)κ(θctrl,iθctrl,i+1)+τjoint,i

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

Iαi(t+Δt)=cωi1(t+Δt)(β+2c)ωi(t+Δt)+cωi+1(t+Δt)+κ(si1(t+Δt)2si(t+Δt)+si+1(t+Δt))κ(θctrl,i(t+Δt)θctrl,i+1(t+Δt))+τjoint,i(t+Δt)=cωi1(t+Δt)(β+2c)ωi(t+Δt)+cωi+1(t+Δt)+κ((si1(t)+ωi1(t+Δt)Δt)2(si(t)+ωi(t+Δt)Δt)+(si+1(t)+ωi+1(t+Δt)Δt))κ(θctrl,i(t+Δt)θctrl,i+1(t+Δt))+τjoint,i(t+Δt)=(c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+κ(si1(t)2si(t)+si+1(t))κ(θctrl,i(t+Δt)θctrl,i+1(t+Δt))+τjoint,i(t+Δt)

The above formula is impossible to integrate because θctrl,i(t+Δt) , θctrl,i+1(t+Δt) , τjoint,i(t+Δt) are unknown at time t.

|κ((θctrl,i(t+Δt)θctrl,i(t))(θctrl,i+1(t+Δt)θctrl,i+1(t)))+(τjoint,i(t+Δt)τjoint,i(t))||(c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+κ(si1(t)2si(t)+si+1(t))κ(θctrl,i(t+Δt)θctrl,i+1(t+Δt))+τjoint,i(t+Δt)|0

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

Iαi(t+Δt)=(c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+κ(si1(t)2si(t)+si+1(t))κ(θctrl,i(t+Δt)θctrl,i+1(t+Δt))+τjoint,i(t+Δt)(c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+κ(si1(t)2si(t)+si+1(t))κ(θctrl,i(t)θctrl,i+1(t))+τjoint,i(t)=(c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+τκ,i(t)+τjoint,i(t)ωi(t+Δt)=ωi(t)+αi(t+Δt)Δt=ωi(t)+ΔtI((c+κΔt)ωi1(t+Δt)(β+2c+2κΔt)ωi(t+Δt)+(c+κΔt)ωi+1(t+Δt)+τκ,i(t)+τjoint,i(t))

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

[ωi1(t+Δt)ωi(t+Δt)ωi+1(t+Δt)][ωi1(t)ωi(t)ωi+1(t)]+ΔtI[(c+κΔt)(β+2c+2κΔt)(c+κΔt)][ωi1(t+Δt)ωi(t+Δt)ωi+1(t+Δt)]+ΔtI[τκ,i1(t)+τjoint,i1(t)τκ,i(t)+τjoint,i(t)τκ,i+1(t)+τjoint,i+1(t)]

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

ω(t)[ωi1(t)ωi(t)ωi+1(t)]
Pn×n[(c+κΔt)(β+2c+2κΔt)(c+κΔt)]
τrem(t)[τκ,i1(t)+τjoint,i1(t)τκ,i(t)+τjoint,i(t)τκ,i+1(t)+τjoint,i+1(t)]

Then, the matrix formula is expressed as follows.

ω(t+Δt)ω(t)+ΔtIPn×nω(t+Δt)+ΔtIτrem(t)

Now, the following approximation can be obtained where In×n is a unit matrix of size n×n.

ω(t+Δt)(In×nΔtIPn×n)1(ω(t)+ΔtIτrem(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×106sec, the calculation error accumulated for the rotational inertia of the whole worm significantly influenced the calculation result ωi and si . To prevent this, the error is corrected as follows. If x¯ixixc , y¯iyiyc , the moment of inertia of the entire worm at time t is as follows by the parallel axis theorem.

Ibody(t)=mi((x¯i(t))2+(y¯i(t))2)+nI

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)

τbody(t+Δt)i(x¯i(t)Fb,i,y(t)y¯i(t)Fb,i,x(t))1+bΔtM

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

Lbody(t)mi(x¯i(t)vi,y(t)y¯i(t)vi,x(t))

The predicted value of ωi at t+Δt, ωip , is calculated by the semi-implicit Euler method (See ‘Numerical integration of the rotational motion’ in Appendix). (ωp=[ωip]T)

ωp(In×nΔtIPn×n)1(ω(t)+ΔtIτrem(t))

The predicted value of si at t+Δt is as follows.

sip=si(t)+ωipΔt

Predicted values xip , yip , vip for xit+Δt , yit+Δt , vit+Δt are calculated from dc , vc , sip , ωip (See ‘Minimum information required to describe the motion of each rod’ in Appendix).

Using xip , yip , vip , the moment of inertia Ibodyp and rotational inertia Lbodyp of the entire worm at time t+Δt are calculated as follows. (where x¯ipxipxcp , y¯ipyipycp , vi=[vi,xvi,y]T)

Ibodyp=mi((x¯ip)2+(y¯ip)2)+nI
Lbodyp=mi(x¯ipvi,ypy¯ipvi,xp)

ωi(t+Δt) is calculated as follows.

ωi(t+Δt)=ωip+(LbodypLbody(t))+τbody(t+Δt)ΔtIbodyp+Ibody(t)2

si(t+Δt) is calculated as follows.

si(t+Δt)=si(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×106sec, as if Δt1×106sec.

Proper selection of friction coefficients

When using the vertical and horizontal friction coefficients bagar, , bagar, 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η,=ηbagar,  , bη,=ηbagar,  by multiplying scaling factor η to the agar friction coefficients bagar,  , bagar,  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 {θctrl,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)θctrl,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 θctrl,i(t) and the response θi(t) . As η decreased from 1 to 10-2 , the total traveled distance of the worm (t|vc(t)Δt|) and the total absolute angle change (S=t=0TΔt|si(t+Δt)isi(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 |Fb,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 {θctrl,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 |Fb,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η,=bagar, /bagar, ), 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
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
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
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)θctrl,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=t=0TΔt|si(t+Δt)isi(t)i|, where T is the total time of the experimental video).

Appendix 1—figure 4
The magnitude of the frictional force |Fb,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 |Fb,i(t)|.

Appendix 1—figure 5
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)θctrl,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=t=0TΔt|si(t+Δt)isi(t)i|)

Appendix 1—figure 6
The magnitude of the frictional force |Fb,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 |Fb,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
    1. Broekmans OD
    2. Rodgers JB
    3. Ryu WS
    4. Stephens GJ
    (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. Book
    1. Hall DH
    2. Altun ZF
    (2008)
    C Elegans Atlas
    Cold Spring Harbor Laboratory Press.
    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
  2. Conference
    1. Lam SK
    2. Pitrou A
    3. Seibert S
    (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
    1. van Aken AA
    2. de Baat C
    3. van Rossum GM
    4. Mulder J
    5. Kalk W
    (1995)
    “Prosthetic condition” and satisfaction with dentures
    Nederlands Tijdschrift Voor Tandheelkunde 102:12–14.
    1. White JG
    2. Southgate E
    3. Thomson JN
    4. Brenner S
    (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

Article and author information

Author details

  1. 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
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0004-2335-3500
  2. 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
  3. 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
    ORCID icon "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

  1. Sent for peer review:
  2. Preprint posted:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. 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

  • 757
    views
  • 75
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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)

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

  1. Taegon Chung
  2. Iksoo Chang
  3. 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

Share this article

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

Further reading

    1. Computational and Systems Biology
    David B Blumenthal, Marta Lucchetta ... Martin H Schaefer
    Research Article

    Degree distributions in protein-protein interaction (PPI) networks are believed to follow a power law (PL). However, technical and study bias affect the experimental procedures for detecting PPIs. For instance, cancer-associated proteins have received disproportional attention. Moreover, bait proteins in large-scale experiments tend to have many false-positive interaction partners. Studying the degree distributions of thousands of PPI networks of controlled provenance, we address the question if PL distributions in observed PPI networks could be explained by these biases alone. Our findings are supported by mathematical models and extensive simulations and indicate that study bias and technical bias suffice to produce the observed PL distribution. It is, hence, problematic to derive hypotheses about the topology of the true biological interactome from the PL distributions in observed PPI networks. Our study casts doubt on the use of the PL property of biological networks as a modeling assumption or quality criterion in network biology.

    1. Computational and Systems Biology
    2. Microbiology and Infectious Disease
    Priya M Christensen, Jonathan Martin ... Kelli L Palmer
    Research Article

    Bacterial membranes are complex and dynamic, arising from an array of evolutionary pressures. One enzyme that alters membrane compositions through covalent lipid modification is MprF. We recently identified that Streptococcus agalactiae MprF synthesizes lysyl-phosphatidylglycerol (Lys-PG) from anionic PG, and a novel cationic lipid, lysyl-glucosyl-diacylglycerol (Lys-Glc-DAG), from neutral glycolipid Glc-DAG. This unexpected result prompted us to investigate whether Lys-Glc-DAG occurs in other MprF-containing bacteria, and whether other novel MprF products exist. Here, we studied protein sequence features determining MprF substrate specificity. First, pairwise analyses identified several streptococcal MprFs synthesizing Lys-Glc-DAG. Second, a restricted Boltzmann machine-guided approach led us to discover an entirely new substrate for MprF in Enterococcus, diglucosyl-diacylglycerol (Glc2-DAG), and an expanded set of organisms that modify glycolipid substrates using MprF. Overall, we combined the wealth of available sequence data with machine learning to model evolutionary constraints on MprF sequences across the bacterial domain, thereby identifying a novel cationic lipid.