Forward and backward locomotion patterns in C. elegans generated by a connectome-based model simulation

Caenorhabditis elegans (C. elegans) can produce various motion patterns despite having only 69 motor neurons and 95 muscle cells. Previous studies successfully elucidate the connectome and role of the respective motor neuron classes related to movement. However, these models have not analyzed the distribution of the synaptic and gap connection weights. In this study, we examined whether a motor neuron and muscle network can generate oscillations for both forward and backward movement and analyzed the distribution of the trained synaptic and gap connection weights through a machine learning approach. This paper presents a connectome-based neural network model consisting of motor neurons of classes A, B, D, AS, and muscle, considering both synaptic and gap connections. A supervised learning method called backpropagation through time was adapted to train the connection parameters by feeding teacher data composed of the command neuron input and muscle cell activation. Simulation results confirmed that the motor neuron circuit could generate oscillations with different phase patterns corresponding to forward and backward movement, and could be switched at arbitrary times according to the binary inputs simulating the output of command neurons. Subsequently, we confirmed that the trained synaptic and gap connection weights followed a Boltzmann-type distribution. It should be noted that the proposed model can be trained to reproduce the activity patterns measured for an animal (HRB4 strain). Therefore, the supervised learning approach adopted in this study may allow further analysis of complex activity patterns associated with movements.

. Diagram of neuro-muscular connections illustrated based on the WormAtlas database 6 . Synaptic connections and gap junctions of whole muscle cells (total 95) and motor neurons VA, VD, VB, DA, DD, and DB. The green oval represents the body-wall muscles, the left side of the figure is the anterior direction, the right side is the posterior direction, and, starting from the upper row, it shows the left dorsal side, right dorsal side, right ventral side, and left ventral side. The circles represent motor neurons, and the difference in color corresponds to the types of motor neurons, as shown in the legend. Solid lines represent synaptic connections, and dotted lines represent gap junctions. www.nature.com/scientificreports/ This study examines whether the motor neuron and muscle network could generate oscillation patterns for both forward and backward movement, and analyzes the distribution of the synapse and gap connection weights. The proposed model consists of motor neurons of classes A, B, D, AS, and muscle. Both gap junctions and synaptic connections are included based on the connectome data. The input to the motor neurons from the command neurons is modeled as a binary input based on their role in evoking forward and backward movements 1,8 , in addition to the generated teacher data from measured muscle activities of a strain expressing calcium indicator in the body-wall muscle (HRB4). The connection parameters were then trained using a supervised algorithm. Finally, the trained model was driven by the untrained input of the command neuron to examine its generalizability, and the fit of trained weights to the modified Boltzmann distribution was tested.

Materials and methods
Assumptions. Figure 1 shows the connection diagram of the model based on the connectome described in the nematode database WormAtlas 6 . In particular, the model consists of seven types of 69 motor neurons that receive five pairs of ten binary inputs, and four rows of 95 muscles 5,6,16 . The neurons and muscles are connected based on the adjacency matrix of the synaptic and gap connections accessible from the following https:// wormw iring. org/ pages/ adjac ency. html 6 .
We considered the inputs to motor neurons from command neurons of PVCL/R and AVBL/R (responsible for forward movement), in addition to AVAL/R, AVDL/R, and AVEL/R (responsible for backward movement) 1 . We assumed that these inputs have no specific polarity, and do not include an oscillating component to generate periodic muscle contractions [17][18][19][20] .
Motor neurons are comprised of DB, DA, DD, VB, VD, VA, and AS (indicated by circles in Fig. 1). Cholinergic neurons, DB, DA, VB, VA, and AS were assumed to excite VD, DD, and muscles. GABAergic neurons VD and DD then inhibit the cholinergic neurons and muscles. Although the polarities of the synaptic connections were determined based on the description in a previous study 5,6 , they have not been fully confirmed experimentally. Therefore, it should be noted that these connections have the potential to be either excitatory or inhibitory in an actual animal.
Proprioceptive feedback also plays an important role in generating undulatory movements, and the undifferentiated processes extending from A-and B-class motor neurons are responsible for proprioceptive feedback 7 . Because the model does not contain the body in order to focus on the neuromuscular system, we connected several muscles to an anterior A-class neuron and a posterior B-class neuron to represent proprioceptive feedback, assuming that muscle activity correlates with body curvature. Here, we did not impose the polarities of the proprioceptive feedback connections.
In Fig. 1, the C. elegans muscle indicated by the ellipse is composed of 95 cells arranged in four rows, and the ventral and dorsal muscle cells alternately contract and relax to produce a smooth crawling motion. For simplification, we assumed the simultaneous activation in each left and right muscle pair, considering their two-dimensional movement on the agar. We trained the model using backpropagation through time (BPTT), a supervised algorithm, to generate muscle activity patterns using two types of teacher data: a mathematically defined sinusoidal pattern, and an activity pattern measured from a fluorescence strain (HBR4).
Mathematical expression of the model. The model considers inputs from forward command neurons (PVCL/R, AVBL/R) and backward command neurons (AVAL/R, AVBL/R, AVEL/R), which are modeled to output control signals of 1 or 0. An output of 1 represents the activated state, and 0 represents the resting state. As shown in the following equations, the forward command neuron L a (t) outputs 1 and the backward command neuron L c (t) outputs 0 during a preset time [T 2v−1 , T 2v ) (v = 1, 2, . . . , o) to command forward movement. This output is reversed during [T 2v , T 2v+1 ) to command backward movement: where L a (t) is the output of the forward command neurons (PVCL/R, AVBL/R) and L c (t) is the output of the backward command neurons (AVAL/R, AVBL/R, AVEL/R). The index a = 1, 2, 3, 4 corresponds to a total of four forward command neurons, and c = 5, 6, …, 10 corresponds to a total of six backward command neurons.
In the motor neurons, a muscle contraction pattern is generated based on the binary input from the command neurons. The motor neuron is defined by the following equations: www.nature.com/scientificreports/ where x N i (t) and x M u (t) are currents flowing into motor neuron i and muscle u at time t, respectively; y N i (t) and y M u (t) are membrane potentials of motor neurons and muscles, respectively; ω l0 i is the bias, τ i is the first-order lag element, and F s is the sampling frequency configured in the simulation procedure. ω NN ij , ω MN iu , ω IN ip are weights representing the synaptic connection strength among motor neurons, that from muscle to motor neuron, and that from the command neuron to motor neuron, respectively.g NN ij = g NN ji , g MN iu = g MN ui are weights corresponding to the conductance of the gap junctions. J = 69 is the number of motor neurons, P = 10 is the number of command neurons, and U = 95 is the number of muscle cells. I p,i (t) represents the proprioceptive feedback to the A-and B-class motor neurons determined by the following equation.
Here, w MN i{3(i+1)−u} denotes the connection strength, 3(i + 1) − u and 3(i + 1) + u are the indexes of a muscle, and we set the proprioceptive feedback length as S = 7 muscles. For the motor neurons located at the head and tail, S was set to the number of the existing anterior and posterior muscles. In this configuration, an A-and B-class neuron receives the outputs from up to seven anterior and posterior muscles.
When the muscles receive the output from the motor neurons, they output the muscle activation level according to the following equations: where ω NM uj is a weight representing the synaptic connection strength from the motor neuron to muscle, and g MM uk is a weight representing the conductance of the gap junctions between the muscle, K = 95 is the number of muscle cells.
According to the connectome data 6 , the absence of connections was represented by restricting the weights of the synapse connection and gap junction to 0 during the following training process.
Training algorithm. The parameters of the neural network model (first-order lag element, synaptic connections, and gap junctions) are trained using the BPTT algorithm 21 . The evaluation function is defined as: is the activity level of the muscles, and T is the maximum step number of the simulation time. BPTT updates each parameter in a particular direction to minimize the evaluation function, which is determined by partial differentiation of the evaluation function for the corresponding parameter. In addition, the following constraint condition is set for each parameter based on the type of connections.
・The synaptic connection strengths from the excitatory neurons DA ( i = 1, 2, . . . ・The first-order lag elements τ i ,τ u are constrained to positive values ( τ i ≥ 0). The teacher data d u (t) used for parameter adjustment are generated based on the motion analysis of C. elegans. As shown in Fig. 2, we tracked seven points on the body of C. elegans while performing forward movement and plotted the time change of the angles between the two straight lines connecting adjacent points. In the figure, each angle shows a sinusoidal-like waveform that is transmitted from the head to the tail. Because the movement of C. elegans is generated by alternating activities of the muscles, the activity level is considered to have similar characteristics. Therefore, we express the teacher data for the activity level of the muscles according to Eqs. (9)- (11).
In the first period of forward movement (T 1 ≤ t < T 2 ) , when the index q = 1, 2, …, 24 is sequentially defined for the muscles from the head to the tail, the teacher data for the left dorsal (u = 1, 2, …, 24) and the right dorsal (u = 25, 26, …, 48) muscles are given by the following equation:  Based on measured data, the angular frequency is set to ω = 1.6 π. Figure 3 shows the procedure of preparing the second teacher dataset. For d u (t) of the second teacher dataset, the muscle activities were measured using a strain expressing calcium indicator in the body-wall muscle (HBR4: goeIs3[pmyo-3::GCamP3.35::unc-54-3' utr, unc-119(+)]V), and the model was trained using the measured fluorescence rate. The specification of the fluorescence microscope system was same as that described in the literature 22 and the locomotion was video-recorded based on procedures 3 previously established for wild-type animals. The recorded video images were then trimmed to extract the duration of one cycle of forward movement, followed by one cycle of backward movement. As shown in Fig. 3, a video analysis software (WormLab, MBF Bioscience, Williston, USA) was used to track the body outline of the animal, and then the body was divided into 24 equal segments from the head to the tail. The fluorescence intensities at the dorsal (k = d) and ventral (k = v) sides of the divided i-th segment F k i,t were measured in the ventral and dorsal areas of each segment, and converted to fluorescence rates R k i is the minimum fluorescence intensity at the k side of the i-th segment. R k i,t was concatenated to generate several cycles of forward and backward movement as shown in Fig. 4, and fitted to multiple sinusoidal functions, as expressed by the following equation, to smooth the measured data: where n = 8 is the number of the functions, t is the time, a ui , ω ui , and φ ui are the amplitude, frequency, and phase of the i-th sinusoidal function, respectively. The fitted data were then concatenated to generate three cycles of (9) d u (t) = sin (ωt + φ u (t)) (u = 1, 2, . . . , 24, 25, 26, . . . , 48), d u (t) = sin (ωt − π + φ u (t)) (u = 49, 50, . . . , 72, 73, 74, . . . , 95).

Simulation conditions
The teacher data for training was configured as follows: the simulation time was set to T = 30 s, and the sampling frequency was set to F s = 0.05 s, the time for switching between forward and backward movement was set to  Figure 4 A illustrates the teacher data generated using Eqs. (1), (2), and (9), and B shows the teacher data generated using the measured muscle activity. The same teacher data were used to train the ipsilateral side, assuming movement on the 2D agar plane. The left part of Fig. 4 shows the binary input from the command neurons. Color maps show the teacher data of the muscles. This figure shows that during forward movement, waves are transmitted from the head to the tail with some phase delay. Conversely, the waves are transmitted from the tail to the head during the period of backward movement. The phase of each wave is reversed at the time when the forward and backward movement switches.

Results
Training of model parameters. Figure 5 illustrates the convergence of the residual errors defined by the evaluation function (Eq. (8)) over ten training trials using the teacher data generated from Eqs. (9)-(11) and another ten training trials using teacher data generated from the measured muscle activity. By adjusting the model parameters, the evaluation function decreased as the learning iterations increased, and finally reached an error tolerance of E ≤ 0.005, although the convergence speeds differ depending on the random initial parameters. Figure 6 shows the residual errors of pre-training, post-training, and verification trials. Significant differences (p < 0.001) between pre-training and post-training, and pre-training and verification trials indicate that the model was successfully trained. Although significant differences were also found between the post-training and verification trials, we confirmed that the model could generate the oscillation of muscle activity and switch between oscillation patterns for forward and backward movements, as described below.
Muscle activity generated by the trained model. Figure 7 shows examples of muscle activities generated by feeding trained input to the model. The two sample data shown in Fig. 7 are those with the median values of residual errors among the ten trials of training, respectively, for the two types of teacher data. The figure confirms that the muscle activity patterns corresponding to forward and backward movement are successfully generated based on the inputs from the command neurons. The muscle activity propagates from head to tail when the forward command neurons activate, and it propagates from tail to head when the backward command neurons activate. However, the model failed to generate activities in neck muscles, as shown in the output of the four-eight muscles from the head in each quadrant. This is because the model focuses on the ventral cord motor neurons and does not include the nerve ring motor neurons. The nerve ring is a neuropile involved in processing various information related to various behaviors in addition to forward and backward movement 5 . Four muscles from the head in each quadrant are innervated by the nerve ring motor neurons, and the next four muscles each in the head are innervated by both the body-wall motor neurons and head motor neurons 23 . This result indicates that the gap junctions between the muscles cannot transmit the activity pattern to the four muscles in the head, and these muscle cells are independently controlled by the nerve ring motor neurons.
In the verification simulation, we inputted command neuron signals to evoke switching between forward and backward movement at different points in time from those in the teacher data. The switching times are T 1 = 0.0, T 2 = 4.7, T 3 = 7.3, T 4 = 17.25, T 5 = 22.5, and T 6 = 28 s for the model trained by the teacher data generated by Eqs. (9)- (11), and T 1 = 0.0, T 2 = 10.1 s for the model trained by the teacher data generated from the measured muscle activity. The outputs of the muscles are shown in Fig. 8. Again, the two sample data shown in Fig. 8 are those with the median values of residual errors among the ten training trials for the two types of teacher data. The   (11). (B) Teacher data generated using the measured muscle activity. The red arrows represent the time duration of forward movement, and the blue arrows represent the time duration of backward movement. The output switches between forward and backward movement at the predetermined time. The x-axis represents the muscle cell number, indexed from head to tail, the y-axis represents time, and the color represents the teacher data for the muscle activation level. The same teacher data are used for the right side of the body-wall muscles. www.nature.com/scientificreports/ results confirm that muscle activity patterns corresponding to forward and backward movement are successfully generated, as with the training data. This result indicates the generalizability of the time of switching between the forward and backward movements.  www.nature.com/scientificreports/ where A and a are the scaling factors, β is the power index, and n is the total number of synaptic connections. We adopted this equation to fit the strength distribution of the synaptic and conductance weights. Figure 9 plots the distribution of the trained synaptic and conductance weights for each of the 10 training trials. The frequencies of weights in all training trials showed a similar trend of decreasing exponentially with the weight strength. www.nature.com/scientificreports/ the modified Boltzmann distribution. In addition, Fig. 9(b) predicts that the conductance of the gap connection could also follow the modified Boltzmann distribution. Theoretically, synapse connections following the modified Boltzmann distribution form a candidate network structures that can maximize the information storage 15 . In this study, the neural network model was trained to store the dynamic patterns corresponding to forward and backward motion using the BPTT algorithm, thus, as a consequence, the weights may have followed a modified Boltzmann distribution.

Discussion
This study formulated a neural network model consisting of 69 motor neurons and 95 muscle cells based on the connectome of C. elegans. After the parameters of the model were adjusted using machine learning, it was shown that the motor neurons could generate the activity patterns of the muscle cells corresponding to forward and backward movement based on the input from the command neurons modeled as binary values. However, the four-eight muscles in the head failed to generate oscillation patterns. This result suggests that the gap junction between muscle cells cannot transmit activity to generate head motion, and the nerve ring motor neuron is necessary to control the activity of these muscles.
Previous models assigned the values of parameters by various means [11][12][13][14][24][25][26][27][28] , such as manually 10 or by using an evolutional algorithm 13,26 . However, the distribution of the assigned parameter values was not analyzed. In this study, we trained the proposed model using the BPTT algorithm. We found that the modified Boltzmann distribution was well fitted to the distribution of the trained parameters (synaptic and conductance weights). This result predicts that the motor neuron and muscle network downstream of the command neurons could form a sparse network in terms of connection strength, which is advantageous to code various movement patterns 15 .
While the previous models reproduced the frequency and wavelength of muscle activities, the proposed model reproduced the fluorescence rates measured from the body wall muscles using a fluorescence strain (HBR4). Therefore, the supervised learning approach taken in this study may allow further analysis of complex activity patterns associated with movements because it provides a framework to reproduce the measured muscle activity patterns of an actual animal.

Conclusion
In this paper, we presented a connectome-based neuromuscular network model of C. elgans. The model is trained to generate muscle oscillation patterns for both backward and forward movements using a supervised learning approach. The main finding of this study is that a motor neuron and muscle network with a sparse connection strength can generate the oscillatory patterns. In addition, the model can be trained to generate measured muscle activity patterns. Therefore, the supervised learning approach taken in this study may allow further analysis of complex activity patterns associated with movements.