A computational model of internal representations of chemical gradients in environments for chemotaxis of Caenorhabditis elegans

The small roundworm Caenorhabditis elegans employs two strategies, termed pirouette and weathervane, which are closely related to the internal representation of chemical gradients parallel and perpendicular to the travelling direction, respectively, to perform chemotaxis. These gradients must be calculated from the chemical information obtained at a single point, because the sensory neurons are located close to each other at the nose tip. To formulate the relationship between this sensory input and internal representations of the chemical gradient, this study proposes a simple computational model derived from the directional decomposition of the chemical concentration at the nose tip that can generate internal representations of the chemical gradient. The ability of the computational model was verified by using a chemotaxis simulator that can simulate the body motions of pirouette and weathervane, which confirmed that the computational model enables the conversion of the sensory input and head-bending angles into both types of gradients with high correlations of approximately r > 0.90 (p < 0.01) with the true gradients. In addition, the chemotaxis index of the model was 0.64, which is slightly higher than that in the actual animal (0.57). In addition, simulation using a connectome-based neural network model confirmed that the proposed computational model is implementable in the actual network structure.

in close proximity at the nose tip. Therefore, the neural network is required to calculate and internally represent two types of gradients using the temporal responses of the sensory neurons. The input to a pair of interneurons (AIYL and AIYR) may play this role, as suggested by a study by Kocabas et al. 6 , which found that the symmetric input to the neurons changes the frequency of pirouette and that the asymmetric input controls the curving rate. In addition, model studies have predicted neural processing related to the weathervane strategy. Ferrée et al. 7 constructed a simple nonlinear neural-network model to enable weathervane, and also extracted its computational rules using the impulse response of a linear neural network 8 . Morse et al. 9 then verified the neural network model by navigating a robot to a light source. Izquierdo and Lockery 10 subsequently proposed a simplified network structure that could adjust the curving rate from sensory inputs, and found that the nonlinearity and self-feedback of motor neurons may serve as key mechanisms for this function. Furthermore, their group demonstrated that the neural network model derived from the connectome was also able to perform weathervane 11 . Xu et al. 12 proposed a dynamic neural-network model and simulated attraction and avoidance behaviour; in addition, their group recently combined a body dynamics model and a neural network model to enable chemotaxis 13 .
An understanding of chemotaxis acquired from the previous studies can be interpreted from the viewpoint of Marr's level of analysis 14 . The behavioural analysis 4,5 highlighted the problem at the computational level that the chemotaxis is closely related to the chemical gradients in the environments. The neural activity measurements 6 revealed the phenomena at the implementational level. The simulation approaches 7-12 analysed the information-processing mechanism at the algorithmic and implementational levels. However, these approaches did not directly treat the algorithm to calculate chemical gradients in the environments, which was explicitly or implicitly assumed given in the computational level. In other words, the relationships between the sensory input at the nose tip of C. elegans and internal representations of the two types of chemical gradients in the environments closely related to the observed strategies have not been formulated explicitly. Therefore, a gap exists between computational and implementational levels. Further, most previous models focused on revealing the mechanism of the weathervane and did not simulate the pirouette simultaneously.
This paper presents a simple and comprehensive computational model based on the motion of the animal involved in chemotaxis to bridge the gap of understanding between computational and implementational levels. The ability of the computational model to convert sensory inputs at the nose tip of C. elegans into internal representations of the chemical gradient parallel and perpendicular to the travelling direction was verified using a chemotaxis simulator that can simulate the body motions of pirouette and weathervane. The chemotaxis performance of the model was compared with previous experimental data 5 . For additional analysis of the implementational level, a connectome-based neural network model was constructed to test whether the computation could be implemented. Based on these results, we also discuss the relationship between the proposed computational model and the findings of the experimental 6 and model 11 studies.

Results
Chemical gradient with respect to the travelling direction derived from the NaCl concentration at the nose tip. Behavioural experiments revealed that pirouette and weathervane strategies are closely related to the chemical gradients parallel and perpendicular to the travelling direction 4,5 , respectively, but how can these gradients related to the traveling direction be obtained by using only the information accessible by the animal? To answer this question, we focused on the fact that ASEL/R responds to the time derivative of NaCl concentration, and the time derivative can be approximated by a directional derivative (see S2 Appendix 2 for the detailed derivation process). The internal representation of the gradients can then be obtained by decomposing the directional derivative of NaCl concentration sensed at the nose tip into the directional components parallel and perpendicular to the travelling direction. Based on this idea, a computational model is derived by using the directional derivative and first mean value theorem for definite integrals. The derived computational model is expressed by equations (1) and (2), to describe the relationships between the chemical concentration sensed at the nose tip and the internal representations of the chemical gradient in the environments, respectively.
where y p and y w are the internal representations of the gradients parallel and perpendicular to the travelling direction at the body centre, respectively, and q 0 is the head-bending angle (cf. Materials and Methods; Multibody model); a p and a w are the reciprocals of the time constants that smooth the time derivative of the NaCl concentration sensed at the nose tip, and b p and b w are the gain constants to scale the inputs. The parameters in the model were adjusted and then set as follows: a p = 0.58, a w = 0.73, b p = 1.20, b w = 1.46. Because the derivation process involves several assumptions and approximations, we use p  and w  to express the accumulated approximation errors (see S2 Appendix 2).
The computational model can be interpreted as follows: Equation (1) works as a low-pass filter to eliminate the head-bending component, and equation (2) allows the comparison of the NaCl concentration between the ventral and dorsal sides using head-bending angles required to calculate the internal representation of the gradient perpendicular to the travelling direction (note that C. elegans lies on its side and navigates by dorsoventral motions). To test the computational model, a chemotaxis simulation involving both weathervane and pirouette strategies was performed by using a multibody model of C. elegans and the chemical environmental model. Chemotaxis simulation using the multibody model. To evaluate the computational model represented by equations (1) and (2), it is necessary to determine the true chemical concentration at the nose tip and responses of sensory neurons (ASEL and ASER) as well as internal representations of the gradient. However, it is difficult to retrieve this information using experimental approaches. Therefore, we constructed a chemotaxis simulator including a multibody model that can perform motions related to pirouette and weathervane and an environmental model that can simulate the NaCl distribution on the agar plate (cf. Materials and methods, Chemotaxis simulator). This method enables the verification of the computational model at the behavioural level without measurement of internal representations of the gradients in the actual neural network. Figure 1a shows travelling path examples of the body centre obtained by chemotaxis simulation using the chemotaxis simulator. The chemotaxis simulation was repeated 10 times. In this simulation, the internal representation of the gradients parallel and perpendicular to the travelling direction were calculated using the computational model (equations (1) and (2)), and the calculated values were used to control the motion of the body model to perform pirouette and weathervane. The average travelling speed of the model was 1.28 ± 0.09 mm/s, which is consistent with the experimental data 5 . Figure 1b compares the chemotaxis index between the body model and the animals. Based on the previously defined chemotaxis index 5 , it is calculated as (T in − T out )/T total , where T in is the time spent within π (2/ ) cm from any NaCl peak, T out is time spent outside the area, and T total is the total simulation and experiment time, which is 1200 s in this case. The figure confirms a chemotaxis index of 0.64 ± 1.38 for the simulation. A comparison between the multibody model and the animal is presented in Supplemental Information S3 from the following aspects: the pirouette motion of the multibody model, a relationship between the curving rate and the gradient perpendicular to the traveling direction, a relationship between the bearing angles and the sharp turn angle, and evaluation of the error caused by converting body postures into the travelling path by using the multibody model.
The true values of the gradients parallel and perpendicular to the travelling direction, denoted as y p G and y w G , respectively, were then geometrically calculated based on the chemotaxis simulation path of the body centre, using equations (14)- (17), as described in Materials and Methods (Geometrical calculations of NaCl gradient), and compared with the internal representation of the gradients calculated using the computational model (equations (1) and (2)). Figure 2  (p < 0.01) and r = 0.89 (p < 0.01) for the gradients parallel and perpendicular to the travelling direction, respectively. The average correlations over 10 chemotaxis simulations were r = 0.90 ± 0.03 and r = 0.91 ± 0.02, and the average root mean square error (RMSE) values were 1.71 ± 0.32 × 10 −6 mM/s and 0.14 ± 0.01 × 10 −3 mM/cm, for the respective gradients. The computational model was derived based on several assumptions and approximations regarding the inputs. We thus tested input (q 0 , dc(x 0 , t)/dt) dependency of the errors, as shown in Fig. 3. Here, the errors are defined as = − e y y p p G p and = − e y y w w G w , respectively, for the gradient parallel and perpendicular to the travelling direction. The figure confirms low correlations between the inputs and the errors. The average correlations of the 10 simulation results are as follows: • q 0 vs e p : 0.01 ± 0.02 • dc(x 0 , t)/dt vs e p : 0.08 ± 0.05 • q 0 vs e w : 0.01 ± 0.02 • dc(x 0 , t)/dt vs e w : 0.1 ± 0.06 *p < 0.001 for all correlations. These results demonstrate low correlations between the errors and inputs, indicating almost no input dependency.  (1) and (2). Spikes observed in the time differential of the NaCl concentration are caused by backward movement for the initiation of pirouette. (b) Comparison of the gradient parallel to the travelling direction. The correlation is r = 0.91 (p < 0.01) and RMSE is 1.78 × 10 −6 mM/s. (c) Comparison of the gradient perpendicular to the travelling direction. The correlation is r = 0.89 (p < 0.01) and RMSE is 0.15 × 10 −3 mM/cm. Generating internal representations of the chemical gradient based on the neural network structure of C. elegans. The computational model indicates that generating the internal representation of NaCl gradients requires the time derivative of NaCl concentration at the nose tip and head-bending angle. We thus tested whether the connectome-based neural network exhibited the ability to generate the internal representations from this information. A neural network model was constructed based on the neural connection structure derived from WormAtlas 15 and neurons included in the model are were chosen based on those treated in Iino and Yoshida 5 . The simulation results were then compared with the computational model described by equations (1) and (2). The structure of the network model is shown in Fig. 4. The model considers the response characteristics of sensory neurons ASEL and ASER using equations (20)-(21) (cf. Materials and Methods, Neural network model) based on measured data 16 . In addition, although previous studies 4, 6,7 suggested that the gradients are not directly coded in the neural network of C. elegans, to facilitate the evaluation of the implementability of the computational model in the neural network, we assumed that the interneurons explicitly represent these gradients. In this case, PVC(L/R) and AVB(L/R) interneurons were assumed to output the gradient parallel to the travelling direction because these neurons are responsible for controlling forward motion 17 and may inhibit backward motion followed by a sharp turn. In addition, the RIA(L/R) and AIZ(L/R) interneurons were assumed to output the gradient parallel to the travelling direction because these neurons control head bending 1,18 , which can generate turning bias. The model considers both chemical synapse connections and gap junctions; these connection weight parameters were trained using back-propagation through a time algorithm 19 modified to simultaneously train the parameters of both gap junctions and chemical synapses. The training datasets were generated using data from the chemotaxis simulation performed using the multibody model, wherein the inputs were the NaCl concentration sensed at the nose tip of the multibody model and the head-bending angles, and the teacher signals comprised true chemical gradients that were geometrically calculated from the path of the body centre obtained from the chemotaxis simulation.  (2)) were calculated for the 10 chemotaxis-simulation datasets using the previously described parameters. Figure 6a,b show examples of the outputs of the neural-network model and the computational model. Figure 6c shows the correlations between the outputs of the neural network model and those of the computational model where the average correlations were r = 0.93 ± 0.04 for both gradients parallel and perpendicular to the travelling direction and p < 0.01 for all simulation datasets.

Discussion
Relationship between NaCl concentration at the nose tip and the internal representations of chemical gradients. This paper presents the computational model that can describe the relationships between the chemical concentration sensed at the nose tip of C. elegans and the internal representations of chemical gradients closely related to pirouette and weathervane strategies, and tests if the computational model is implementable in the connectome-based neural network model. This analysis concept is based on Marr's level of To verify the computational model, we constructed a chemotaxis simulator which includes the environment model expressing NaCl distribution and the multibody model of C. elegans. The results shown in Fig. 1 and the supplemental data indicate that the chemotaxis simulator has the ability to simulate the behaviour during chemotaxis, and is therefore applicable for testing equations (1) and (2). Figure 2 then confirms that equations (1) and (2) could convert the NaCl concentration sensed at the nose tip into the internal representation of the NaCl gradient parallel and perpendicular to the travelling direction. Equation (2) indicates that generating the internal representation of the gradient perpendicular to the travelling direction requires the head-bending angle. It is commonly considered that head bending is generated by a central pattern generator 20 involving stretch receptors 21 ; thus, its information may be obtained through these neurons.
The derivation process (see S2 Appendix 2) shows that the directional decomposition of the time derivative of NaCl concentration requires cosine (symmetric) and sine (asymmetric) functions of the head-bending angle, respectively, for parallel and perpendicular to the travelling direction. Similar decomposition process could be performed in the pair of AIY neurons or their upstream neurons because the results of an experimental study 6 indicated that the symmetric input of the pair of AIY neurons controls the pirouette frequency and the asymmetric input controls the gradual turn (weathervane).
In addition, one of the findings of the model study by Izquierdo & Beer 11 indicated that the asymmetrical response characteristics of the ventral and dorsal motor neurons are required to perform gradual turns (weathervane). That is, the response characteristic of either side of the motor neuron would be shifted to the region of lower sensitivity and the other side to that of higher sensitivity, by the sensory input, so that the motor neurons generate a biased sinusoidal wave to regulate body motion. From the computational model (equation (2)), it can be interpreted that the motor neurons perform the directional decomposition, and the gradient parallel to the travelling direction is coded in the motion of the animal.
As described above, the proposed computational model suggests that both findings obtained from the observation of AIY neurons and motor neurons in the simulation could be explained by directional decomposition. From the viewpoint of Marr's level of analysis, the computational model can bridge the gap between the problem defined in the computational level and implementational level.
Generating internal representations of chemical gradients based on the neural network structure of C. elegans. Figure 5 confirms that the training algorithm successfully adjusted the parameters of  the neural network model and that it could internally represent gradients with respect to the travelling direction using the NaCl concentration sensed at the nose tip and the head-bending angle. High correlations between the respective gradients generated by the computational model, neural network model, and geometrical calculation (shown in Fig. 6) indicate that the computational model can be implemented in the connectome-based neural network.

Limitation.
In the neural network model, AVB and PVC were assumed to output the gradient parallel to the travelling direction. However, in the measurement experiments, activities in these neurons were observed despite the absence of a gradient 22 , which indicates that AVB and PVC are redundant in the calculation of the gradients. We confirmed this redundancy by excluding AVB and PVC from the connectome-based neural network model, setting AIY as the output neuron, and readjusting the parameters. Accordingly, AIY generated the gradient parallel to the travelling direction with a correlation of 0.98 (p < 0.001), indicating that the computational model could be implemented with various parameter sets. This redundancy in the parameter space limits the analysis of the contribution of the respective neuron to chemotaxis. The parameters of the multibody model and motions such as body forms and reversal durations are chosen such that the body model can approximate the motions of the actual animal as close as possible (see Supplemental Information S3). Although the approximation error affects the analysis results, the error of generated gradients (see Fig. 2) caused by the multibody model may be partial (see Fig. S5 in Supplemental Information S3). The modelling error, however, is inevitable in a process to approximate the motion of C. elegans. To reduce the modelling error, merging the simulation and experimental approaches can be a candidate solution for future research. In addition, the motion of the multibody model was given by a sinusoidal function, as shown in equations (5) and (6) in Material and Methods: Chemotaxis simulator. Thus, the model ignores small random wriggling motions of the nose tip that can be observed in the actual animal. As the sensory neurons are located at the nose tip, this motion can affect the sensory inputs. However, this evaluation is not performed in this study.

Conclusion
In this study, a simple and comprehensive computational model was derived to convert the response of a single sensory input into two types of internal representations of the NaCl gradient parallel and perpendicular to the travelling direction and enabled simultaneous simulation of the pirouette and weathervane strategies. The derived computational model suggests that internal representations of the gradients can be generated by combining head-bending angles and sensory input from ASEL/R neurons. It could also be used to interpret the functions of AIY neurons and motor neurons, respectively, identified in previous experimental 6 and simulation studies 11 , and thus can bridge the gap between the chemotaxis problems at the computational and implementational levels.
The connectome-based neural network model included in the chemotaxis simulator demonstrated that the computational model could be implemented in it, although the coding manner of the chemical gradient might differ from that of the actual animals. The connectome-based neural network model may allow further analysis of the functions of respective neurons by introducing the biological constrictions and measured neural activities and by simulating ablation experiments.

Materials and Methods
To evaluate the derived computational model [equations (1) and (2)], we developed a chemotaxis simulator including the multibody and environmental models. The body model can perform both weathervane and pirouette motion. The NaCl chemical gradients parallel and perpendicular to the travelling direction in an environment, which are called a true gradient in this paper, were calculated from the simulated body centre path. To test the implementability of the proposed computational model under actual neural structure representations, a connectome-based neural network model was constructed. All simulations described in this study were performed using MATLAB 2013. Model parameters are given in S1 Appendix 1. The following text explains each part of the chemotaxis simulator, simulation method, and neural network model.
Chemotaxis simulator: Multibody model. C. elegans body was approximated using a multibody model, as shown in Fig. 7, defined by the following Newton-Euler equations based on a previous study 23 : are centrifugal and Coriolis force terms, g(q) and g g (q) are gravity force terms, ρ = [ρ 0 , ρ 1 , …, ρ L−1 ] T is the driving torque generated from the motors, F j = [F T,j , F N,j ] is the friction vector between the floor and the j-th body, and J j is the Jacobian matrix. Please note that we omit the explicit notation of time dependence of the variables for simplification, but all variables depend on time in equations (3) and (4). The parameters related to inertia, centrifugal, Coriolis, and gravity forces are determined based on animal size and weight (Tables in S1 Appendix 1). The friction forces were determined based on the average velocity of the animal on the chemotaxis plate (0.12 mm/s) 5 .
Chemotaxis simulator: implementing pirouette and weathervane. Body model motion is controlled by a sinusoidal function determining each i-th joint angle q i (t).

Amphid and ASE neurons
i where q Max is the maximum joint bending angle, ω determines the angular velocity of bending, φ i (t) determines the wavelength of the animal body as well as switching between forward and backward motion, and ψ (in rad) is determined by the measured body wavelength 24 . The bias parameter κ(t) = κ w (t) + κ R (t) is determined by the weathervane turning rate κ w (t) and random walk κ R (t). Pirouette and weathervane were implemented by changing φ i (t) and κ(t) as follows: Weathervane.
w ws where d s (t) (in mM/mm) represents the gradient perpendicular to the traveling direction and C w is derived from the measured weathervane curvature 5 .
Pirouette. Pirouette frequency f p (t) (in events/s) is given by the following equation 5 : where d t (t) (in mM/s) is the gradient parallel to the travelling direction. A uniform random number R(t) ∈ [0, 1] is generated every sampling time and pirouette is initiated if R(t) < f p (t).
The pirouette starts with backward motion by switching the phase φ i (t) using the following equation: where T p0 is the pirouette initiation time. The backward motion lasts for T b = 6.0 s because measurement data indicate that most sharp turns occur after a long reversal with around three head-bending cycles 25 and their probability is approximately 90% after 6.0 s of backward motion 26 . When the backward motion is finished, equation (9) is applied again to switch back to forward motion. A sharp turn is then initiated after the backward motion using the following equation: where T Ω0 = T p0 + T b is the sharp turn initiation time, T Ω1 = T Ω0 + T p1 , T Ω2 = T Ω1 + T p2 , and T ΩE = T Ω2 + T p3 are duration parameters for linearly changing the φ i (t), and C p is the maximum phase delay. The total sharp turn duration is set to T ΩE − T Ω0 = 3.18 s, and C p was adjusted to generate turns larger than 1.75 rad based on definitions from a previous study 5 . Temporal changes in ϕ i (t) allow the body shape to change from the S shape to the Ω shape and back, as illustrated in Fig. S2 in Supplemental Information S3.
Random walk. The random walk was introduced based on measured curvature in the absence of a chemical gradient as in the following equations: T r Based on simulations performed by Iino and Yoshida 5 , a random curvature parameter b T was chosen every ΔT = 12 s and linearly altered for every sampling time Δt (in s).
Chemotaxis simulator: Environmental model. The environmental model and its parameters were derived from a previous study 5 . NaCl concentration at an arbitrary position x can be calculated by solving Fick's equation as follows: where N 0 is the NaCl solution concentration, D is the NaCl diffusion coefficient, E is the agar plate thickness, and x k is the coordinate of the k-th NaCl point. where x g (t) − x g (t − Δt) is the body centre travelling direction and δ is a small displacement; thus Δx v (t) and Δx d (t) are the vectors perpendicular to the travelling direction.
Neural network model. The neural network model shown in Fig. 4 was defined based on the actual connection structure derived from WormAtlas 15 considering both synaptic connections and gap junction, and neurons included in the model were derived from Iino and Yoshida 5 . The model receives inputs of the head-bending angles and NaCl concentration at the nose tip of the multibody model from RIV(L/R) and amphid sensory neurons ASE(L/R), respectively. The ASEL and ASER responses are determined from the previous experimental data 16 and sent to the interneurons in both amphid and nerve ring. To facilitate the evaluation of implementability of the computational model on the neural network, we assumed that the interneurons could directly generate the gradients, although previous studies 6, 10,11 have suggested that these gradients are not directly coded in the C. elegans neural network. In this case, PVC(L/R) and AVB(L/R) interneurons were assumed to output the gradient parallel to the travelling direction because these neurons are responsible for controlling forward motion 17 and may inhibit backward motion followed by a sharp turn. Additionally, interneurons RIA(L/R) and AIZ(L/R) were assumed to output the gradient perpendicular to the travelling direction because these neurons control head bending 1,13 , which can generate turning bias. The following equations give a mathematical definition of the neural network model: The dynamic characteristic of a neuron is given as follows: where u i (t) corresponds to the electrical current to the i-th neuron, U i (t) corresponds to the membrane potential, T s is the sampling time, τ i is the time constant caused by leakage current, w ij , and w ij d are the chemical synapse connection strengths from the j-th to the i-th neuron, s ij = s ji is the conductance of a gap junction, n is the total number of neurons, and m is the total input number. Additionally, z j (t) = u j (t) when j ≤ n, and z j (t) = I j−n (t) when n < j ≤ n + m where I j (t) is the external input including NaCl concentration and head-bending angle q 0 (t).
ASER and ASEL neurons were modelled based on the experimental data 16  where I 1 (t + 1) and I 2 (t + 1) represent inputs to ASEL and ASER neurons, respectively. dc/dt is the NaCl concentration time derivative at the nose tip. The parameters R R R R are adjusted to fit the response peaks of ASEL and ASER neurons. Figure 8 shows the fitting results, which confirms that the model could generate a response similar to the experimental data 16 .
Adjustments of chemical synapse connections w ij and gap junctions s ij were performed using back-propagation through a time algorithm, which uses the chain-rule of partial differential on the following evaluation function: where μ i = 1 if index i corresponds to the neuron that outputs a gradient either parallel or perpendicular to the traveling direction; otherwise μ i = 0, d i (t) is the target gradient that is desired to be outputted from the corresponding neuron. The chemical gradient is normalised to produce d i (t) by using the following equations as where c ζ is a gradient either parallel or perpendicular to the traveling direction. The parameters w ij , s ij , τ i are then iteratively updated based on the partial differential of H by the following equations: