MetaWorm: An Integrative Data-Driven Model Simulating <i>C. elegans</i> Brain, Body and Environment Interactions — part 5/6

Cleaned-text code points [59157, 73491). Bibliography omitted here; source.md includes it.

Our neural network model consists of 136 neurons, including 15 sensory neurons for inputs 
and 80 motor neurons for outputs (Supplement Table S5). Here are the reasons why we chose 
these neurons. First, 65 head ganglia neurons, whose neural dynamics were recorded in the 
published whole- brain Ca2+ imaging experiment 17 should be selected to enable the 
comparison between experiment and simulation. Second, other 71 neurons, most of which 
were motor neurons in the locomotion circuit of C. elegans7,18, were added to the network so 
that the model could be complete in a sensory-motor functional structure. In summary, this 
model contained neurons for chemosensory, decision, and locomotion. 
1.5 Connection parameters optimization 
The weights of synapses and gap junctions in the neural network was initialized randomly. 
The polarities of synapses were initialized randomly according to a given 
excitatory/inhibitory ratio (2.33). These parameters were optimized using an unpublished 
optimization algorithm for multi-compartment neural network. The algorithm employed a 
mathematical approximation approach to calculate gradients for Back Propagation Through 
Time (BPTT) to update parameters. Except for the weights (of synapses and gap junctions) 
and the polarities (of synapses), the inputs (of some sensory neurons) also can be optimized 
by this optimization algorithm. 
1.6 PCA of neural activities 
Principal Component Analysis (PCA) was applied to the time series data of soma membrane 
potential of the 65 head neurons. Before conducting the analysis, the time series data was 
preprocessed by subtracting the mean value of the trace. 
1.7 Correlation analysis  
The correlation between two neurons was quantified by the Pearson correlation coefficient 
of their neural activities, which ranged from -1 to 1. A coefficient of -1 or 1 indicates a strong 
negative or positive correlation of neural activity, respectively. A coefficient of 0 indicates 
that the neural activities of the two neurons are unrelated. In this analysis, the neural activity 
was represented by the membrane potential of soma. 
2 Details of Body & Environment Model 
2.1 The body model 
Tetrahedron is the simplest ordinary convex polyhedra in three-dimensional space. It has 
four vertices, six edges, and is bounded by four triangular faces. We used tetrahedron as the 
basic element to construct the body model. Surface shape of the body was extracted as 
triangle list from Sibernetic’ membrane data19. Then we triangulated the body with a robust 
and high-quality algorithm20. The triangulated tetrahedron mesh was refined for latter 
mechanics solver. The final body mesh was in a straight, natural resting pose that we defined 
as the standard body (Fig.6A). The triangulated tetrahedron body was not strictly symmetric 
in dorsal-ventral direction and left-right direction. Since triangulation algorithm had 
additional constraints on quality of tetrahedrons. This non-strict symmetry is probably more 
consistent with real C. elegans. 
2.2 Muscles 
96
 muscles were arranged from head to tail according to C. elegans  anatomy. Each 24 
muscles were grouped into longitudinal string, marked as ventral right (VR), ventral left 
(VL), dorsal right (DR), dorsal left (DL) in four quadrants (Fig.6B)19. We imported body 
membrane mesh into blender21 and modeled VR, VL, DR, DL by four Bezier curves. Then 
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 
The copyright holder for this preprintthis version posted March 10, 2024. ; https://doi.org/10.1101/2024.02.22.581686doi: bioRxiv preprint

muscle geometry was exported as wavefront obj file. We used Finite Element Method (FEM) 
to simulate the soft body, and used finite element constraints to simulate elastic muscle cells 
of body. Each muscle cell (Fig.6A) corresponds to one or more FEM constraints according to 
the discretized tetrahedron body (Fig.6C). Only nearby tetrahedrons can be driven by 
corresponding muscle cells.  
2.3 Soft body solver  
We used projective dynamics22 as the solver of the soft body. The purpose of mechanics 
solver is calculating positions and velocities of all 3D vertices of the deformed body over 
time. Key idea of projective dynamics is to compute a local step and a global step iteratively. 
At local step, vertex positions of the body are fixed and projection variables are updated. At 
global step, projections are fixed and vertex positions of the body are updated. The local step 
is well suited for parallelism, thus can be optimized for performance. The global step is a 
linear system and can be precomputed if the body geometry and constraints are not changed. 
After a few local-global iterations the vertex positions of deformed body are converged. Then 
velocity of deformed body can be computed by Euler integration.  
2.4 Interaction with environment 
A simplified model of hydrodynamics was used to simulate the interaction between the 
body and the fluid, similar to the published study 23. For simplicity, we supposed external 
force is only formed at the surface of body. Two kinds of external forces were considered 
here: thrust force offered by body movement, and drag force offered by fluid. With this 
simplification, computation in the whole 3D fluid environment is avoided. Thus, more agents, 
larger environments, more complex tasks are achievable through this simplification. 
2.5 Numerical stable coordinate system for body movement 
C. elegans  is a soft body animal which has special patterns during movement. Every 
position of the body is oscillating during movement. Moreover, there is no apparent important 
marks on the body during locomotion, therefore it’s hard to find an apparent stable coordinate 
system for C. elegans movement. Simply using head or body center as a coordinate system 
can’t analysis the body’s movement quantitatively and stably. Here we proposed an effective 
and simple method to measure C. elegans locomotion. The core of our method is to find a 
numerically stable coordinate system relative to the body itself. The main steps of our method 
are as follows:  
1) Define the target body, standard body and standard body coordinate system (SBCS).  
2) Obtain transformation between the target body and the standard body. 
3) Use the transformation in 2) to perform three-dimensional transformation on the SBCS to 
obtain the target body reference coordinate system (TBRCS). 
Details of the three steps are described below: 
Step 1) The target body is the research subjects. It can change its shape over time and 
exhibits various behaviors. The standard body is a special case of the target body. Here we 
define the standard body in a natural state where the head and tail are in a straight line, 
without stretching and bending. There is no internal force in the standard body soft body, and 
no external force acting on the standard body surface. The S BCS is a three-dimensional 
orthogonal coordinate system used to reflect the overall structure and function of the standard 
body. In this case, we define the SBCS with origin at the standard body center. X-axis of the 
SBCS is towards the head of the standard body. Z- axis of the S BCS is towards ventral to 
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 
The copyright holder for this preprintthis version posted March 10, 2024. ; https://doi.org/10.1101/2024.02.22.581686doi: bioRxiv preprint

dorsal direction of the standard body. Three axes of the SBCS form a right-handed coordinate 
system. Schematic diagram of the target body, the standard body and the SBCS can been seen 
in Figure. 7A. 
Step 2) In this step, we first define the corresponding matching points on the target body 
and the standard body. For example, simply use all the vertices of the body. The 
transformation between the matching points on the target body and the standard body is 
calculated as the transformation between the target body and the standard body. Calculation 
formulas of the transformation are as follows: 
𝑝𝑝0��� = 1
𝑚𝑚 � 𝑝𝑝𝑖𝑖
0
𝑚𝑚
𝑖𝑖=1
 
𝑝𝑝̅ = 1
𝑚𝑚 � 𝑝𝑝𝑖𝑖
𝑚𝑚
𝑖𝑖=1
 
𝑞𝑞𝑖𝑖
0 = 𝑝𝑝𝑖𝑖
0 − 𝑝𝑝0��� 
𝑞𝑞𝑖𝑖 = 𝑝𝑝𝑖𝑖 − 𝑝𝑝̅ 
𝐻𝐻 = 𝑞𝑞𝑖𝑖
0𝑞𝑞𝑖𝑖𝑇𝑇 
𝑅𝑅 = 𝑉𝑉𝑈𝑈𝑇𝑇 
𝑇𝑇 = 𝑝𝑝̅ − 𝑅𝑅𝑝𝑝 0��� 
𝑀𝑀 = �𝑅𝑅 𝑇𝑇
0 1 � 
The matching points are represented by column vectors. 𝑝𝑝𝑖𝑖
0  are matching points on the 
standard body, and 𝑝𝑝𝑖𝑖 are matching points on the target body . m is the number of matching 
points. 𝑝𝑝0��� is the center point of matching points on the standard body , and 𝑝𝑝̅ is the center 
point of matching points on the standard body. 𝑞𝑞𝑖𝑖
0 refers to the standard body matching points 
with the center offset subtracted, and 𝑞𝑞𝑖𝑖 denotes the standard body matching points with the 
center offset subtracted. 𝐻𝐻  is a matrix computed by 𝑞𝑞𝑖𝑖
0𝑞𝑞𝑖𝑖𝑇𝑇 . 𝑞𝑞𝑖𝑖𝑇𝑇  is transpose of 𝑞𝑞𝑖𝑖 , 𝑞𝑞𝑖𝑖 =
[𝑞𝑞1, 𝑞𝑞2 … 𝑞𝑞𝑚𝑚], 𝑞𝑞𝑖𝑖
0 = [𝑞𝑞1
0, 𝑞𝑞2
0 … 𝑞𝑞𝑚𝑚0 ]. 𝑈𝑈 and 𝑉𝑉  are obtained by the singular value decomposition 
𝐻𝐻 = 𝑈𝑈Λ𝑉𝑉 of the matrix  𝐻𝐻. 𝑈𝑈𝑇𝑇 is transpose of 𝑈𝑈. 𝑀𝑀 is the transformation matrix between the 
target body and the standard body, which will be used in next step. 𝑅𝑅 is rotation part of 𝑀𝑀, 𝑇𝑇 
is translation part of 𝑀𝑀. A schematic diagram of matching points and transformation 𝑀𝑀 are 
shown in Fig. 7A. 
Step 3) In this step, we use matrix 𝑀𝑀 to perform three -dimensional transformation on the 
SBCS to obtain the TBRCS . The matrix of the TBRCS can be formulated as: 𝑀𝑀� = 𝑀𝑀𝑀𝑀0 =
�𝑅𝑅𝑅𝑅0 𝑅𝑅𝑇𝑇0 + 𝑇𝑇
0 1 � . 𝑀𝑀 = �𝑅𝑅 𝑇𝑇
0 1 � ,  𝑅𝑅  and 𝑇𝑇  have been described in last section. 𝑀𝑀0 =
�𝑅𝑅0 𝑇𝑇0
0 1 � is matrix used to represent the S BCS. 𝑇𝑇0  is the origin of the S BCS, and 𝑅𝑅0 
indicates orientation of the S BCS. The column vector of 𝑅𝑅0  represents the orthogonal 
coordinate axis. Finally, we can make the TBRCS through 𝑀𝑀�: 𝑅𝑅𝑇𝑇0 + 𝑇𝑇 as the origin and 𝑅𝑅𝑅𝑅0 
as the orthogonal coordinate axis.  
2.6 Numerically stable measurements of body movement 
1) Measurement of body trajectory 
During behavior analysis, C. elegans  can be seen as a deformable biology composed of 
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 
The copyright holder for this preprintthis version posted March 10, 2024. ; https://doi.org/10.1101/2024.02.22.581686doi: bioRxiv preprint

multiple points. When studying the locomotion of a deformable body, generally one or more 
points of biology body is defined as tracking object, and locomotion of deformable body is 
studied through tracking these defined points. For deformable biology like C. elegans , 
position and velocity of local tracking points are hard to reflect overall state of body 
locomotion (Fig. 7C). As a result, values such as body  trajectory or steering angle obtained 
through these tracking points are usually unstable over time. Using the TBRCS, we can better 
describe the overall trend of body movement. Comparison of trajectories are shown in Figure. 
7A. If the center point of body is used as trajectory, the trajectory is a zigzag shape. However, 
if the TBRCS is used as trajectory of body movement, the trajectory is a smooth curve which 
means stable values. It can be seen the TBRCS  is more suitable as a measurement of overall 
path and direction of body locomotion. 
2) Tracking body movement 
H
ere, we demonstrate how to track body movement with the TBRCS. If the standard body 
is drawn in the TBRCS , then it can be used as a reference for the target body . The standard 
body and the target body will be center alignment, translate and rotate together. We define 𝑤𝑤𝑘𝑘
0 
as 3D tracking points of the standard body , and 𝑤𝑤𝑘𝑘 as 3D tracking points of the target body. 
k = 1, 2 … n. n is the number of tracking point. Then relative movement of body can be 
defined as 𝑤𝑤𝑘𝑘 − 𝑤𝑤𝑘𝑘
0. It’s a measurement of body locomotion relative to itself (the standard 
body), which is more numerically stable than relative to world. 𝑤𝑤𝑘𝑘 − 𝑤𝑤𝑘𝑘
0 is also better than 
relative to body center, since both translation part and rotation part of whole deformation 
body are considered in the TBRCS. Moreover, sampling bias of tracking points is eliminated 
through subtract operation, since bias minus bias equals zero. 
In
 Figure 7B, we sampled 17 tracking points on the body uniformly. Offset of the tracking 
points are indicated by pink lines between tracking points. Relative velocity of 3D tracking 
points is drawn in blue arrows. As a common method, trajectories of sampled 17 tracking 
points that relate to world are also shown in Fig ure 7C. We can see that compared to all 17 
tracking points trajectories, trajectory of the TBRCS  is the smoothest. Relative position and 
relative velocity of 3D tracking points over time are plotted in Figure 7 E and Fig ure 7 F 
respectively. In Figure 7E, Positions of 3D tracking points are tiled from bottom to top for 
clear seeing. In this case, swing direction (Z direction of the TBRCS ) of swimming body can 
be clearly reflected with the TBRCS . We can also find that there is a spring movement in X 
direction of the TBRCS (tail to head) during forward swimming. Locomotion in Y direction 
of the TBRCS is nearly flat since in this case muscle activations are approximately 
symmetrical in this direction. From Fig ure 7 F we can further deduced that the biggest 
propulsion of forward swimming is come from tail (purple line), since speed means power. 
3) Steering angle of body movement 
If 𝑂𝑂𝑡𝑡𝑖𝑖   is the trajectory of TBRCS origin at time 𝑑𝑑𝑖𝑖 . 𝑑𝑑𝑖𝑖 =  𝑑𝑑1, 𝑑𝑑2 … 𝑑𝑑𝑚𝑚 . Then we can get 
velocity of TBRCS  by 𝑣𝑣𝑖𝑖 =
𝑂𝑂𝑡𝑡𝑖𝑖−𝑂𝑂𝑡𝑡𝑖𝑖−1
𝑡𝑡𝑖𝑖−𝑡𝑡𝑖𝑖−1
. In this way, velocity magnitude and velocity direction 
over time are numerically stable. We used angle θ which is between TBRCS x-axis and 𝑣𝑣𝑖𝑖 as 
a measurement of steering, as shown in Fig ure 7D. Both 𝑣𝑣𝑖𝑖 and the TBRCS  are numerically 
stable over time, thus θ is also numerically stable. 
3 Locomotion Behavior Simulation of C. elegans  
3.1 Simulation configuration 
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 
The copyright holder for this preprintthis version posted March 10, 2024. ; https://doi.org/10.1101/2024.02.22.581686doi: bioRxiv preprint