A neural network model that generates salt concentration memory-dependent chemotaxis in Caenorhabditis elegans

Two types of extensive evolutionary searches for the electrophysiological parameters of the neuroanatomical minimal network model were performed by optimizing a chemotaxis index (CI). The first approach was a conventional evolutionary search without additional neurophysiological constraints, as had been employed previously (Izquierdo and Beer, 2013). The second evolutionary search was conducted with the constraint that the AIY–AIZ synaptic connections be inhibitory (Li et al., 2014). In these genetic algorithms, the networks evolved in chemotaxis assays where conical shapes were employed as salt concentration profiles (Equation A15). In contrast, the CI for optimized networks was evaluated using a Gaussian shape of salt concentration profile (Equation A16), which mimics the salt concentration gradients observed in laboratory tests of chemotaxis in C. elegans (Ferrée and Lockery, 1999; Ward, 1973). The CI values for the most optimized networks with and without the neurophysiological constraints were 0.877 and 0.855, respectively, as shown in Figure 2. These values were found to be comparable to those reported in the previous study (Izquierdo and Beer, 2013). The difference between these CI values is slight, while the model optimized with the constraints exhibits a marginally accelerated attainment of the salt concentration peak, as shown by the trajectories. The slightly higher chemotaxis performance observed in the constrained model is not essentially attributed to the introduction of the AIY–AIZ synaptic constraints but rather depends on the specific individuals selected from the optimized individuals obtained from the evolutionary algorithm. In fact, even when the AIY–AIZ constraints are taken into consideration, the model retains a significant degree of freedom to reproduce salt klinotaxis due to the presence of a substantial parameter space. Consequently, the impact of the AIY–AIZ constraints on the optimization of the CI is expected to be negligible.

Figure 2

Download asset Open asset

The locomotion trajectories of the most optimized network with the constraints are depicted for the cases in which the model worm was oriented at ten different angles at the initial position (Figure 2b; see also Video 1). For comparison, the trajectories of the most optimized network without the constraints are also depicted in Figure 2a. These trajectories show that both the model worms rapidly turned in the direction of the steepest gradient. These model worms exhibited the salt preference behavior in klinotaxis that were consistent with the experimental observations (Iino and Yoshida, 2009) and previous simulation studies (Chen et al., 2022; Izquierdo and Beer, 2013; Izquierdo and Lockery, 2010).

Video 1

Download asset

Subsequently, in order to determine how the signal arising from alterations in salt concentration was transmitted through the neural circuit to regulate motor systems, we examined the impact of step changes in salt concentration of varying magnitude with different sign on the neurotransmitter release $z_i$ (defined by Equation A5) from each neuron. Figure 3a, b presents the sign and strength of the weight of the synaptic connection, ${\displaystyle w_{ij}}$, in the unconstrained and constrained most optimized models employed in this analysis, respectively. Figure 3c, d illustrates the results of the analysis of $z_i$ for the most optimized networks where the additional neurophysiological constraints are not considered and are considered, respectively. It is noted that $z_i$ differs from the membrane potential $y_i$ as a consequence of the influence of the bias term ${\displaystyle {\theta }_i}$ in Equation A5. Moreover, the neuronal activities involved in signal transmission across chemical synapses are primarily detected through $z_i$. In contrast, $y_i$ correlates with the neuronal activity that are involved in gap junction connections, which can be observed experimentally through an increase in intracellular Ca²⁺ concentration. As illustrated in Figure 3b, the evolutionary search with the constraints successfully yielded the network with inhibitory AIY–AIZ synaptic connections as the most optimized model. Coincidentally, the evolutionary search without the constraints also yielded the network with inhibitory AIY–AIZ synaptic connections as the most optimized model (Figure 3a). However, this is an incidental finding, as the unconstrained evolutionary search also yielded high-CI networks with excitatory AIY–AIZ synaptic connections, which were analogous to the network that had been presented as a representative model presented in the previous study (Izquierdo and Beer, 2013). A noteworthy distinction between Figure 3a and b is that the ASER–AIY synaptic connections in the constrained model are inhibitory (Figure 3b) whereas those in the unconstrained model are excitatory (Figure 3a).

Figure 3 with 3 supplements see all

Download asset Open asset

Related to these, the signal transmissions between the neurons that are involved in the neural circuit for klinotaxis have been investigated using calcium imaging techniques (Matsumoto et al., 2024). As the salt concentration decreased, the activity of the salt-sensing neuron ASER increased, resulting in an observed activation pattern of the interneuron AIZ that was synchronized with the activation of ASER. As illustrated by the red patterns of $z_i$ in Figure 3d, the constrained model exhibited a depolarization of AIZ interneurons as the signal output from ASER increased during a decrease in the salt concentration. This indicates that the model reproduces the pattern of neuronal activity observed in the experiments. In contrast, the unconstrained model exhibited an inverse activity pattern of AIZs, namely, a hyperpolarization of AIZs, in response to a decrease in the salt concentration (Figure 3c). Nevertheless, the CI performance in the unconstrained model was as high as that in the constrained model. It is therefore essential to consider not only the connectome but also the available neurophysiological information in order to identify potential neural circuits that are consistent with the actual worm among the high-CI networks obtained by the evolutionary search. The inhibitory connections between ASER and AIYs are essential for synchronizing the activation of AIZs with the activity of ASER, provided that the AIY–AIZ connections are inhibitory, as has been identified through experiments. (Li et al., 2014; Matsumoto et al., 2024). The biochemical mechanism and neurophysiological significance of the ASER–AIY inhibitory connections will be discussed in detail in a subsequent section.

Next, it was investigated how the most optimized model with the constraints reaches the peak of the salt gradient during klinotaxis. To elucidate the mechanism of klinotaxis, the curving rate was examined in relation to the bearing (Figure 4a). The curving rate, as defined in Figure 4b, was introduced to quantify the behavior of worms during klinotaxis (Iino and Yoshida, 2009). The bearing is defined in Figure 9 and described in detail in Bearing in Methods. The positive correlation between the curving rate and the bearing, as illustrated in Figure 4a, indicates that the model worm continuously turns in the direction of the steepest increase in salt concentration (see Figure 9). It is, next, necessary to reveal how the model worm determines the direction of its turn in order to increase the salt concentration as it moves. The positive correlation between the curving rate and the normal gradient of the salt concentration, as illustrated in Figure 4c, indicates that the model worm is able to sense the concentration gradient normal to the travel direction and thus turn in the direction that increases the salt concentration. These results are consistent with the mechanism that will be proposed later, which indicates that the model worms efficiently sense a slight concentration gradient normal to the direction of travel based on a change in salt concentration with sweeping of the head sensory neurons due to the undulating motions and utilize this information to turn in the correct direction. Figure 4d illustrates the positive and negative components of curving rate, as well as the average of these components, as a function of the translational gradient of salt concentration. The translational gradient is defined in Normal and translational salt concentration gradient in Methods. The positive and negative components of the curving rate are, respectively, sampled from the trajectory during leftward turns (as illustrated in Figure 4b) and rightward turns, respectively. The average of the positive and negative components of the curving rate was close to zero regardless of the translational gradient, due to the symmetry between the dorsal and ventral regions of the model worm. In contrast, the magnitude of the positive and negative components of the curving rate increased as the translational gradient varied from positive to negative. This indicates that, when the translational gradient is negatively large, the model worm regulates the direction of movement with a large positive or negative curving rate based on the normal gradient. Conversely, when the translational gradient is positive, the worm fine-tunes the direction of movement with a small positive or negative curving rate, which is determined by the normal gradient.

Figure 4

Download asset Open asset

Next, focusing on the most optimized network with the constraints (Li et al., 2014; Matsumoto et al., 2024), the signaling pathways in response to step changes in salt concentration are discussed in detail. Changes in the $z_{\text{ON}}$ signal from the ASEL are transmitted to the AIYL via the excitatory connection (Figure 3d). The signals indicating increases in $z_i$ from the AIYL are transmitted to the AIZL via the inhibitory connection between AIYL and AIZL displayed in Figure 3b, resulting in a hyperpolarization of the membrane potential $y_i$ in the AIZL (see Figure 3—figure supplement 1f). The hyperpolarization signals in the AIZL are transmitted to the AIZR via the gap junction (Figure 3—figure supplement 1d, f, Figure 3d). This is because the neuron dynamics via gap junctions results from the equilibration of the membrane potential y_i of two neurons connected by gap junctions rather than the z_i. The signals indicating decreases in $z_i$ from the AIZR are transmitted to the SMBDR and SMBVR via the excitatory connection (Figure 3d). Note that the oscillatory components $z_i^{\mathrm{O}\mathrm{S}\mathrm{C}}$, as defined by Equations A10a and A10b, are introduced into the SMBD and SMBV, respectively. Consequently, the baseline components of $z_i$ from the SMBDR and SMBVR are shifted by half a cycle from each other (see the right side of the fourth and third panels from the bottom in Figure 3d). The discrepancies between the $z_i$ values from the SMBD and SMBV neurons (see Equation A13a) result in an increase in the turning angle $\phi$ relative to its ideal oscillatory component (black line), as illustrated on the left side of the second panel from the bottom in Figure 3d. The implications of the changes in $\phi$ will be discussed later.

Changes in the $z_{\mathrm{O}\mathrm{F}\mathrm{F}}$ signal from the ASER are transmitted to the AIYL via the inhibitory connection (Figure 3d). The signals indicating decreases in $z_i$ from the AIYL are transmitted to the AIZL via the inhibitory connection displayed in Figure 3b, resulting in a depolarization of the membrane potential $y_i$ in the AIZL (Figure 3d). The signals of the depolarization in the AIZL are transferred to the AIZR via the gap junction (see Figure 3—figure supplement 1d, f, Figure 3d). The signals indicating increases in $z_i$ from the AIZR are transmitted to the SMBDR and SMBVR via the excitatory connection (Figure 3d). The disparities between the $z_i$ values from the SMBD and SMBV neurons (see Equation A13a) result in a decrease in the turning angle $\phi$ relative to its ideal oscillatory component (black line), as illustrated on the right side of the second panel from the bottom in Figure 3d. It is remarkable that although the signaling pathway in the most optimized network without the constraints (Figure 3, Figure 3—figure supplement 1) is largely different from that with the constraints, the regulatory mechanisms of $\phi$ exhibited by these networks are analogous. This is evident from a comparison of the second panels from the bottom in Figure 3c, d. As illustrated by the analysis of these signaling pathways, the transmission of signals via electrical gap junctions plays a crucial role in both the neural circuits. Indeed, it was observed in those neural circuits that blocking the gap junctions led to significant changes in neurotransmitter release $z_i$ and resulted in changes in the turning angle ${\displaystyle \phi }$ (see Figure 3—figure supplement 2).

It was then analyzed how changes in $\phi$ in response to step changes in salt concentration caused the worm to move in the correct direction. From the sinusoidal trajectory of the model worm (the left side of the first panel from the bottom in Figure 3d), it can be observed that the head sensory neurons sweep toward the positive y-axis direction during step increases in ${\displaystyle z_{\text{ON}}}$ occurred. This corresponds to a situation where the salt concentrations on the positive side of the y-axis are greater than those on the negative side. The turning angle $\phi$ is increased from its ideal oscillatory component to a value close to zero, causing the model worm to deviate from the ideal sinusoidal trajectory and gradually turn toward higher salt concentrations. Conversely, in the case of step increases in ${\displaystyle z_{\mathrm{O}\mathrm{F}\mathrm{F}}}$, the head sensory neurons sweep toward the positive y-axis direction, as illustrated in the right lowest panel of Figure 3d. This corresponds to a situation where the salt concentrations on the positive side of the y-axis are lower than those on the negative side. The turning angle $\phi$ is reduced from its ideal oscillatory component to be a more negative value, allowing the model worm to turn faster than the ideal sinusoidal trajectory and thus reach higher salt concentrations with greater efficiency. As illustrated by the second panel from the bottom in Figure 3c, the most optimized network without the constraints also uses a similar regulatory mechanism of $\phi$ to efficiently reach the peak of the salt gradient. To verify the applicability of the elucidated regulatory mechanism of $\phi$ to other situations, the changes in $\phi$ in response to step increases in $z_{\text{ON}}$ and $z_{\text{OFF}}$ that were delayed by half of the cycle were examined as a further representative case. This ensured that the regulatory mechanism was operating as intended (see Figure 3—figure supplement 3). These observations demonstrate that the change in $\phi$ relative to its ideal oscillatory component is properly regulated based on the direction of the sensory neuron sweeps at the time of receiving ON and OFF signals via the sensory neurons. The regulatory mechanism of $\phi$ was found to be essentially identical to that reported in the previous studies (Chen et al., 2022; Izquierdo and Beer, 2013; Izquierdo and Lockery, 2010). It is remarkable that this regulatory mechanism derived via the optimization of the CI has been observed in the context of chemotaxis in Drosophila larvae chemotaxis (Wystrach et al., 2016) and phototaxis in zebrafish (Wolf et al., 2017). The principle of operation, in which the dependence of motor responses to sensory inputs on the phase of motor oscillation, therefore, may serve as a convergent solution for taxis and navigation across species.

The neural networks generated by the evolutionary algorithm using the CI (Equation A17) as the fitness function exhibited klinotaxis, with the model worm turning to move toward higher salt concentrations. This salt preference behavior is consistent with the chemotaxis observed in well-fed individuals cultivated at a higher salt concentration than the current environmental concentrations (C_cult > C_test). Conversely, experimental observations have shown that well-fed individuals that are cultivated at a lower salt concentration than the current environment (C_cult < C_test) exhibit the opposite salt preference behavior (Kunitomo et al., 2013). However, the neural circuit mechanism underlying the reversal of preference behavior in klinotaxis remains unclear.

The reversal of chemotaxis, which depends on the cultivated environmental salt concentration, has been investigated in detail with respect to klinokinesis (Hiroki et al., 2022; Sato et al., 2021). The interneuron postsynaptic to the ASER, AIB, which is involved in the behavior of klinokinesis, exhibits the bidirectional responses to changes in salt concentration in the current environment, depending on the cultivated salt concentration (Sato et al., 2021). The bidirectional responses of the AIB are found to be attributed to a change in the basal level of glutamate neurotransmitter released from the ASER (Sato et al., 2021). Specifically, the differential sensitivities of the excitatory glutamate receptor GLR-1 and the inhibitory glutamate receptor AVR-14, which are expressed on the AIB, result in the bidirectional response of the postsynaptic AIB depending on the changes in the basal level of glutamate release from the ASER (Hiroki et al., 2022). The GLR-1 is a glutamate-gated cation channel, whereas the AVR-14 is a glutamate-gated chloride anion channel. In addition, the sensitivity of GLR-1 (EC₅₀ 5 mM) is several tens of times lower than that of the AVR-14 (EC₅₀ 0.2 mM) (Hiroki et al., 2022). Consequently, the synaptic transmission between the ASER and AIB is excitatory and inhibitory, respectively, when the basal level of glutamate release from the ASER is high and low. In addition, the biochemical mechanism of the change in the basal level of glutamate neurotransmitter released from the ASER was also elucidated. When the cultivated salt concentration exceeds the current environmental concentrations (C_cult > C_test), the basal glutamate release from the ASER is enhanced by the phosphorylation of UNC-64/Syntaxin 1A at Ser65 through the protein kinase C (PKC-1) signaling pathway (Hiroki et al., 2022). Conversely, in the case of C_cult < C_test, the basal glutamate release is decreased as a consequence of the reduced PKC-1 activity (Hiroki et al., 2022). The AIB and the neurons further downstream of the AIB show the salt concentration memory-dependent bidirectional responses that correlate with the reversal of the salt preference behavior in klinokinesis (pirouettes). These observations suggest that the reversal of salt preference behavior in klinokinesis is attributed to the synaptic plasticity, which varies from the excitatory to inhibitory synaptic connections between ASER and AIB.

Given those experimental findings, we propose a potential mechanism that may underlie the salt concentration memory-dependent reversal of preferential behavior in klinotaxis. At least two glutamate receptors, the inhibitory glutamate-gated chloride anion channel GLC-3 (Ohnishi et al., 2011) and the excitatory metabotropic glutamate receptor MGL-1 (Kang and Avery, 2009), have been demonstrated to be expressed in the AIY, the postsynaptic interneuron to the ASER. The EC₅₀ for GLC-3 has been reported to be 1.9 mM (Horoszok et al., 2001), while the EC₅₀ for MGL-1 is not available in the literature. On the other hand, the EC₅₀ for the C. elegans homolog of MGL-1, MGL-2, has been reported to be 9 µM (Tharmalingam et al., 2012). The sensitivity of MGL-2 is ~200 times higher than that of GLC-3, suggesting that the sensitivity of MGL-1 may also be sufficiently higher than that of GLC-3. In the case of C_cult > C_test, if the basal level of glutamate release from the ASER is nearly comparable to the level of EC₅₀ for the GLC-3, the activity of MGL-1 is saturated so that the MGL-1 does not contribute to a change in the membrane potential of the AIY upon the depolarization of the ASER (for more details regarding the release of internal Ca²⁺ from the endoplasmic reticulum, see, e.g., the reference on a metabotropic acetylcholine receptor; Sumi and Harada, 2023). Thus, an inhibitory synaptic transmission between the ASER and the AIY results from the influx of Cl⁻ into the cytosol of the AIY via the GLC-3 upon depolarization of the ASER due to a decrease in salt concentration. Conversely, in the case of C_cult < C_test, when the basal level of glutamate release from the ASER is sufficiently lower than the level of EC₅₀ for the GLC-3 but still nearly comparable to the level of EC₅₀ for the MGL-1, the synaptic transmissions between the ASER and the AIY become excitatory as a result of the Ca²⁺ release from the endoplasmic reticulum, which is mediated by the activation of MGL-1 (Ashida et al., 2019; Shidara et al., 2013). A similar increase in Ca²⁺ in the AIY was also observed in a previous study on a thermotaxis in C. elegans (Ohnishi et al., 2011; Ohta and Kuhara, 2013), although the precise mechanism remains unclear. On the other hand, it is noteworthy that the biochemical mechanism proposed here does not contradict the ASER–AIY inhibitory connections, as identified in the most optimized network with the constraints.

In light of the above argument, we investigated the impact of the salt concentration memory-dependent synaptic plasticity between ASER and AIY on the salt preference behavior in klinotaxis. To this end, we used the most optimized network with the constraints, in which we varied the synaptic connections between ASER and AIY from inhibitory to excitatory. Figure 5a, b illustrates the network model with the inhibitory (equivalent to Figure 3b) and excitatory connections, respectively. The intermediate states between Figure 5a and b are illustrated in Figure 5—figure supplement 1. Figure 5c shows the curving rates obtained from the most optimized network with the inhibitory (Figure 5a) and excitatory (Figure 5b) connections as a function of the normal gradient of salt concentration. For comparison, the curving rate obtained from the ninth intermediate connection (see #9 in Figure 5—figure supplement 1), which is derived by introducing an increment of 1.5 in the synaptic weight ${\displaystyle w_{ji}}$ between ASER and AIY at each step, is also presented in Figure 5c (also shown as #9 in Figure 5—figure supplement 2). Note that the most optimized model with the inhibitory connections (see Figure 5a, also shown as #0 in Figure 5—figure supplement 1) is considered as an individual that is well fed and cultivated at a higher salt concentration than the current environment. The curving rate was varied from an increasing trend function with increasing the normal concentration gradient to a function with a decreasing trend (Figure 5c, see also Figure 5—figure supplement 2) when the weight of the ASER–AIY synaptic connection $w_{ji}$ was increased from a negative (inhibitory) to a positive (excitatory) value. These findings are consistent with the experimental observation shown in Figure 5d, indicating that the most optimized network with the constraints, of which only the inhibitory connections between ASER and AIY were replaced by excitatory connections (Figure 5b), can reproduce the salt preference behavior of a well-fed individual cultivated at a lower salt concentration than the current environment.

Figure 5 with 5 supplements see all

Download asset Open asset

A similar analysis to that displayed in Figure 3d is presented in Figure 5e, with the exception that the inhibitory connections between ASER and AIY have been replaced by the excitatory connections as illustrated in Figure 5b. Note that the changes in $z_i$ and $\phi$ in response to step changes in $z_{\text{ON}}$ (depicted in blue) observed in the neurons further downstream of the ASEL are fully identical to those illustrated in Figure 3d, as the synaptic connections between ASEL and AIY remain unaltered. The direction of changes in $z_i$ and $\phi$ in response to step changes in $z_{\text{OFF}}$, as indicated in red, accords with that of $z_i$ and $\phi$ in response to step changes in $z_{\text{ON}}$, as shown in blue (Figure 5e). This observation indicates that even when either $z_{\text{OFF}}$ or $z_{\text{ON}}$ increases during a head sweep, the $\phi$ is varied from its ideal oscillatory component to zero. This causes the model worm diverging from the ideal sinusoidal trajectory and gradually turning in the direction of the head sweep that occurred during the ON or OFF signal. More specifically, the model worm gradually turns in the direction of either the lower or higher side of the salt concentration upon receiving the OFF or ON signal, respectively (see both the first and second panels from the bottom in Figure 5, Figure 5—figure supplement 3). Nevertheless, as evidenced by the curving rate in Figure 5c, the model worm clearly shows a preference for lower salt concentrations in klinotaxis. Therefore, the behavior of turning to lower salt concentrations upon increasing $z_{\text{OFF}}$ should play a more decisive role in the reversal of the salt preference behavior than the behavior of turning toward higher salt concentrations upon increasing $z_{\text{ON}}$. Indeed, as illustrated in Figure 5—figure supplement 4b, the majority of trajectories yielded by the neural circuit shown in Figure 5b demonstrated that the worm model turned to move in the opposite direction of the peak of the salt concentration gradient (see also Video 2). However, a subset exhibited a slight curve toward the gradient peak, which then proceeded straight and passed without any discernible response (Figure 5—figure supplement 4b and Video 2). These behaviors can be interpreted in terms of the changes in the turning angle φ in response to step increases in $z_{\text{ON}}$ and $z_{\text{OFF}}$ (Figure 5, Figure 5—figure supplement 3). In the case of step increases in $z_{\text{OFF}}$ as illustrated in the second right panel from the bottom in Figure 5e, the turning angle φ is increased from its ideal oscillatory component to a value close to zero, causing the model worm to deviate from the ideal sinusoidal trajectory and gradually turn toward lower salt concentrations. On the other hand, in the case of step increases in $z_{\text{ON}}$ as illustrated in the second left panel from the bottom in Figure 5e, the turning angle φ is again increased from its ideal oscillatory component to a value close to zero, causing the model worm to deviate from the ideal sinusoidal trajectory and gradually turn toward higher salt concentrations. The behaviors that are consistent with these analyses are observed in the trajectory illustrated in Figure 5—figure supplement 4b. For a detailed analysis of signal transmission through chemical synaptic connections and electrical gap junctions related to these behaviors, Figure 5—figure supplement 5 is presented. In summary, the reversal of salt preference behavior in klinotaxis of a well-fed individual cultivated at a lower salt concentration than the current concentrations can be primarily attributed to the change in synaptic connections between ASER and AIY from inhibitory to excitatory connections.

Video 2

Download asset

For C. elegans, finding food and dwelling at the food source are essential survival strategies. In the absence of food, the release of dopamine from PDE neurons is suppressed, the AVK interneurons, which make inhibitory connections with the PDE, in turn results in the release of FLP-1 neuropeptides (Oranth et al., 2018). As a result, the FLP-1 neuropeptides inhibit SMB motor neurons (Oranth et al., 2018), which, in conjunction with the reversal of starvation-induced klinokinesis (Kunitomo et al., 2013), may result in the dispersal behavior observed in starved individuals. The dispersal behavior of starved individuals, which is necessary for searching distant food, has been reported to manifest as alterations in salt preference behavior in klinotaxis as well as klinokinesis (Kunitomo et al., 2013). In particular, the salt concentration memory-dependent preference behavior observed in well-fed individuals (see Figure 5d) is markedly suppressed in starved individuals, as illustrated in Figure 6a (Kunitomo et al., 2013). This implies that in order to efficiently reach distant food sources, the worm suppresses turning and tends to move in a straight line. Given the above neurophysiological findings, we investigated the impact of inhibiting SMB activity on the salt preference behavior in klinotaxis. To this end, we employed the neural circuits in which the ASER–AIY transmission was inhibitory, excitatory, and intermediate between the two, that is, the three circuits that yielded the results shown in Figure 5c. Figure 6b shows that the inhibition of the SMB activity by reducing the strength of synaptic connections between AIZ and SMB resulted in the suppression of salt preference behavior, which is consistent with the experimental observations presented in Figure 6a. Moreover, the inhibition of SMB function through the reduction of bias terms in the SMB motor neurons resulted in the suppression of the salt preference behavior (see Figure 6—figure supplement 1b), which is a comparable result to that shown in Figure 6b and is also consistent with the experimental data presented in Figure 6a or Figure 6—figure supplement 1a; Kunitomo et al., 2013. A comparison of Figure 6c, d (or Figure 6—figure supplement 1c, d) with Figure 5—figure supplement 4a, b, respectively, showed that the neural circuit models, in which the SMB functions were suppressed, yielded weak spreading behaviors (see also Video 3).

Figure 6 with 1 supplement see all

Download asset Open asset

Video 3

Download asset

The modeling of ASE sensory neurons does not include either the all-or-none depolarization characteristic of the ASEL or the hyperpolarization characteristic of the ASER in response to step increases in salt concentration (Suzuki et al., 2008). As a result, a comprehensive understanding of the specific roles played by each the ASEL and ASER could not be achieved. In the present study, the oscillator components of the SMB are not intrinsically generated by an oscillator circuit but are instead externally imposed via $z_i^{\mathrm{O}\mathrm{S}\mathrm{C}}$. Furthermore, the complex and context-dependent responses of ASER (Luo et al., 2014) were not taken into consideration. It should be acknowledged as a limitation of this study that these omitted factors may interact with circuit dynamics in ways that are not captured by the current simplified implementation. The available data on synaptic transmission between the AIZ and the SMB, as well as the available information on the activation pattern of the SMB in response to changes in salt concentration, were insufficient to perform the evolutionary search of parameters involving the SMB with the application of appropriate constraints. Consequently, a comprehensive investigation of the role of the SMB in the regulation of the turning direction was hindered by these limitations.

Synaptic transmissions from the chemosensory neurons to interneurons were modeled for the ON (ASEL) and OFF (ASER) sensory neurons, respectively, as an instantaneous function of a coarse-grained time derivative of the recent history of NaCl concentration, $z$ (Izquierdo and Beer, 2013):

where ${\displaystyle c\left(t\right)}$ is the salt concentration at time t, and N and M are the durations of the first and second half intervals, respectively, over which the concentration is averaged. The factor of 100 in Equation A1 was introduced to convert z to values that were more readily amenable to analysis. The ASEL is the ON cell, where the strength of neurotransmitter release $z_{\mathrm{O}\mathrm{N}}$ increases in response to an increase in $c\left(t\right)$. Conversely, the ASER is the OFF cell, with $z_{\mathrm{O}\mathrm{F}\mathrm{F}}$ increasing in response to a decrease in $c\left(t\right)$ (Suzuki et al., 2008; see Figure 7).

Figure 7

Download asset Open asset

The modeling of interneurons used current-based synaptic transmissions as follows Izquierdo and Beer, 2013:

where $y_i$ represents the membrane potential relative to the resting potential of zero, ${\tau }_i$ is the time constant, the first sum represents the input current from chemical synapses, and the second sum represents the input current from electrical synapses. In the first sum, the weight $w_{ji}$ represents the strength of the synaptic connection, and $z_j$ is considered as the strength of neurotransmitter release, as represented by

with the following sigmoid function:

${\displaystyle {\theta }_i}$ in Equation A5 represents the bias term that shifts the range of sensitivity associated with neurotransmitter release. In the second sum of Equation A3, ${\displaystyle g_{ki}}$ represents the conductance of the electrical gap junction, with ${\displaystyle g_{ki}>0}$. The fourth term $I_i^{\mathrm{I}\mathrm{N}}$ in Equation A3 represents the input current from the chemosensory neurons. This is given by

where ${\displaystyle w_i^{\mathrm{O}\mathrm{N}}}$ and ${\displaystyle w_i^{\mathrm{O}\mathrm{F}\mathrm{F}}}$ are the strengths of the synaptic connections between interneuron i and the ASEL and ASER sensory neurons, respectively.

Neck motor neurons were modeled in a manner analogous to that employed for interneurons, except for the inclusion of self-connections (Izquierdo and Beer, 2013):

where the second term represents the input current from chemical synapses, as described in Equation A5. The third term represents a constitutive input current from an oscillatory component,

where $w_{\text{OSC}}$ represents the strength of the connection with the oscillatory component and ${\displaystyle z_i^{\mathrm{O}\mathrm{S}\mathrm{C}}}$ represents the strength of neurotransmitter release from the oscillatory component. Since the dorsal motor and ventral motor neurons receive a constitutive out-of-phase oscillatory input from the motor systems, respectively, it was postulated that modeling ${\displaystyle z_i^{\mathrm{O}\mathrm{S}\mathrm{C}}}$ would be achieved by using a sinusoidal function for the dorsal and ventral motor neurons, respectively:

where $T$ represents the duration of a one cycle of sinusoidal motion on agar, which was estimated to be 4.2 s (Ferrée and Lockery, 1999).

To reduce the total number of parameters included in the model, in accordance with the neuroanatomical minimal network model (Izquierdo and Beer, 2013), the following constraints were introduced with respect to the second term on the right-hand side in Equation A8 for the dorsal and ventral motor neurons on the left and right sides, respectively:

These parameters are constrained to be symmetric between the dorsal and ventral motor neurons on the left and right sides, respectively. As discussed in the previous studies (Izquierdo and Beer, 2013; Izquierdo and Lockery, 2010), a neuron is considered unistable if a self-connection is presented and the weight ${\displaystyle w_i^{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}$ is less than 4. This allows the motor neuron with such a self-connection ${\displaystyle \left(w_i^{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}<4\right)}$ to respond smoothly to different changes in salt concentration. Similarly, the above symmetry constraints were also applied to the weight ${\displaystyle w_i^{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}}$ for the self-connections:

The worm model was constructed using the neuroanatomical minimal network model as the basis for modeling worm locomotion (Izquierdo and Beer, 2013). The locomotion of the worm was modeled as the movement of a single point (r_x, r_y) at a constant velocity v (Figure 8a). The angle of the direction of motion μ relative to the x-axis was measured with counterclockwise being positive (Figure 8b). The model involves two assumptions regarding the locomotion mechanism. The first is that the length of the neck muscles is proportional to the strength of neurotransmitter release from the neck motor neurons (Figure 8a). The second is that the turning angle ${\displaystyle \phi }$ (Figure 8b) is proportional to the difference in the length of the neck muscles. According to these assumptions, the equation of motion governing the angle of movement direction μ is given by

Figure 8

Download asset Open asset

The summations in Equation A13a are performed over the indices ${\displaystyle i_{\text{D}}\in \left\{\mathrm{S}\mathrm{M}\mathrm{B}\mathrm{D}\mathrm{L},\mathrm{S}\mathrm{M}\mathrm{B}\mathrm{D}\mathrm{R}\right\}}$ and $i_{\text{V}}\in \left\{\mathrm{S}\mathrm{M}\mathrm{B}\mathrm{V}\mathrm{L}\text{, SMBVR}\right\}$, respectively. The $w_{NMJ}$ represents the strength of the connection from motor neurons to muscles. In Equation 13b, Equation 13c, $y_{i_{\mathrm{D}}}$ and $y_{i_{\mathrm{V}}}$ represent the membrane potential of the dorsal and ventral neck motor neurons, respectively, and ${\theta }_{i_{\mathrm{D}}}$ and ${\theta }_{i_{\mathrm{V}}}$ represent the bias terms that shift the sensitivity of the neurotransmitter release in the dorsal and ventral neck motor neurons, respectively. $z_{i_{\text{D}}}$ and $z_{i_{\text{V}}}$ represent the strength of neurotransmitter release from the dorsal and ventral neck motor neurons, respectively.

The position (r_x, r_y) of the model worm is updated by the following equation:

where $v$ is a constant velocity of 0.022 cm/s (Ferrée and Lockery, 1999). The present model does not explicitly consider the kinematic mechanism that is responsible for generating the forward thrust due to the oscillatory motion of the worm. However, it is a matter of experimental observation that the locomotion of real worms never occurs without the thrust generated by undulations (Gray and Lissmann, 1964). In order to impose this constraint on the motion of the worm model, the fitness of the individuals that do not exhibit undulations was reduced, thus excluding such individuals from the search performed by the evolutionary algorithm (see also Evaluation of fitness).

Equations A4, A8, A13a, and A14 were solved using the Euler integration with a time step of 0.01 s for the evolution of the network parameters and the evaluation of the chemotaxis behavior for evolved individuals. On the other hand, a time step of 0.001 s was used in the simulations that displayed the trajectories.

The salt concentrations observed during a typical salt chemotaxis assay have a Gaussian distribution (Ward, 1973). However, in the context of evolutionary algorithms, however, Gaussian gradients are problematic because the local steepness depends on the distance from the gradient peak. To circumvent this problem, conical gradients with randomly varying steepness were implemented in each simulation throughout the evolutionary process. Salt concentrations are proportional to the Euclidean distance from the gradient peak ${\displaystyle \left(x_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}},y_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}\right)}$,

where $\alpha$ is the slope of the gradient. Conversely, a Gaussian distribution was used to evaluate the chemotaxis behavior of most optimized individuals:

where $C_0$ = 1.0 (mM) and $\lambda$ = 1.61 (cm) were employed (Ferrée and Lockery, 1999; Ward, 1973).

A genetic algorithm (Back, 1996) was used to evolve the following parameters of the model (Izquierdo and Beer, 2013) in the ranges given in brackets: $w_{NMJ}$ in Equation A13a [1, 3]; all the biases $\theta$ [−15, 15]; the weights $w_{ji}$, $w_i^{\mathrm{O}\mathrm{N}}$ in Equation A7, $w_i^{\mathrm{O}\mathrm{F}\mathrm{F}}$ in Equation A7, and $w_i^{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{f}}$ in Equation A12a [−15, 15]; $w_{OSC}$ in Equation A9 [0, 15], $N$ and $M$ in Equation A1 [0.1, 4.2]; the weight of gap junction $g_{ki}$ in Equation A4 [0, 2.5] (Olivares et al., 2018). In the present study, the genetic algorithm was executed 100 times to obtain 100 optimal networks. The following is a description of a single iteration of the genetic algorithm:

1. An initial population of 60 individuals was randomly generated in which the parameters of the model were encoded in a 22-element vector of real values between [–1, 1], and these parameters were linearly mapped to their corresponding ranges.

2. The fitness of 60 individuals was quantified by using the CI that will be defined by Equation A17. Here, the CI for each individual was obtained as the mean CI value averaged over 50 chemotaxis assays explained later. The top 20 individuals were selected based on the fitness obtained.

3. A next generation population of 60 individuals was generated by three operations:

   1. The selected top 20 individuals were copied as the individuals comprising the next generation population.

   2. The two-point crossover was applied, where the pairs of an odd and the next even number individual in the CI order were selected from the top 20 individuals with a probability of 0.6 as the first and second parents.

   3. Mutations were assumed to occur in the individuals that were selected from the top 20 individuals with a probability of 0.5. The Gaussian noise with a standard deviation of 0.05 was added to the elements selected with a probability of 0.4 from the 22 elements in the vector of the mutants.

   4. If the total number of individuals generated by the selection, crossover, and mutations was less than 60, the remaining individuals were randomly generated in the same way as in procedure 1.

4. Procedures 2 and 3 were repeated to evolve the population for 300 generations.

5. The individual with the highest-CI value after 300 generations of evolution was selected as the best performing individual.

The 100 iterations of this genetic algorithm were performed to generate the ensemble of the 100 best performing individuals. The generated top 100 individuals were ranked based on the CI value.

The fitness of the individuals generated during the evolutionary algorithm was evaluated by performing the simulations with the chemotaxis assays. At the beginning of each simulation, the model worm was placed at the origin ${\displaystyle (x,y)=(0.0,0.0)}$ with an initial orientation randomized over the range ${\displaystyle [0.0,2\pi ]}$ and motor neuron activations were randomized over the range ${\displaystyle [0,1]}$. The gradient steepness $\alpha$ in Equation A15 was randomized over the range ${\displaystyle [-0.38,-0.01]}$ (Izquierdo and Beer, 2013). Fitness was quantified as the CI, defined as the time average of the distance to the peak of the gradient at ${\displaystyle \left(x_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}},y_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}\right)=\left(4.5,0.0\right)}$,

where ${\displaystyle h\left(t\right)}$ is the Euclidean distance from the peak,

$h\left(0\right)$ is the initial distance of the worm from the salt peak, namely, 4.5 cm, and $T_{\mathrm{s}\mathrm{i}\mathrm{m}}$ is the total simulated assay time. $T_{\mathrm{s}\mathrm{i}\mathrm{m}}$ was set to 500 s in the evolutionary algorithm and 1000 s in the evaluation of the CI value for the individuals obtained. For simplicity, negative CI values were set to zero. The fitness of an individual generated during the evolutionary algorithm was determined as the averaged CI values over 50 trials, in which the gradient steepness $\alpha$ of the salt concentration in Equation A15 was randomly chosen in the range of [−0.38,–0.01] and the orientation of the worm at the starting point (0.0, 0.0) was randomized (Izquierdo and Beer, 2013). To ensure the survival of individuals moving in an undulating manner, the CI value was reduced when sinusoidal movements were not present. Specifically, the sign of the turning angle $\phi$ was checked at 1/4 and 3/4 cycles within each period, and if they had the same sign, the CI value was reduced by 0.008 each time as a penalty.

In the calculation of the quantities that were used to characterize the klinotaxis behavior, that is, curving rate, bearing, and salt concentration gradients, the coordinates of the trajectory for the first three cycles of the worm’s sinusoidal locomotion were not used because the time derivative of the NaCl concentration history in Equation A1 could not be precisely determined around the starting point, thus affecting the locomotion. The ensemble average of these quantities was determined by performing 100,000 sets of the simulation with a simulation time of ${\displaystyle T_{\mathrm{s}\mathrm{i}\mathrm{m}}=200\mathrm{s}}$. The Gaussian distribution of the salt concentration given by Equation A16 was used in the analysis of the salt preference behavior in klinotaxis.

Curving rate is also called turning bias (Izquierdo and Beer, 2013). The coordinates at three points over six cycles, namely, zero, $2n\pi$, and $4n\pi$, where $n$ is three, were used to calculate the vectors from the start point to the midpoint $\text{r}\left(2n\pi \right)$ (see Figure 9) and from the midpoint to the end point ${\displaystyle \text{r}\left(4n\pi \right)}$. The turning angle over the six cycles $\phi \left(4n\pi \right)$ was calculated using $\text{r}\left(2n\pi \right)$ and $\text{r}\left(4n\pi \right)$, with counterclockwise being positive (see also Figure 8b, although the time interval between the vectors is different). The locomotion distance $d\left(4n\pi \right)$ was estimated from the sum of $\left|\mathbf{r}\left(2n\pi \right)\right|$ and $\left|\mathbf{r}\left(4n\pi \right)\right|$. The curving rate was determined by ${\displaystyle \phi \left(4n\pi \right)/d\left(4n\pi \right)}$.

Figure 9

Download asset Open asset

The bearing was determined as the angular difference between the vector of translational movement defined by $\text{r}\left(2n\pi \right)$ introduced above and the vector from the current position to the peak of the salt concentration gradient with counterclockwise being positive (see Figure 9).

The normal gradient of salt concentration was calculated from the difference in salt concentration between the position shifted by 0.001 cm toward the left perpendicular to the direction of the worm movement and the current position. The translational gradient of salt concentration was calculated from the difference in salt concentration between the position shifted by 0.001 cm toward the direction of the worm movement and the current position.

Full details of the methods to perform the simulations and the input data have been provided in the main text. The newly developed codes for the simulations and data analysis have been made publicly available on GitHub (copy archived at Hironaka and Sumi, 2025).

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

1. Sent for peer review: November 4, 2024
2. Preprint posted: November 5, 2024
3. Reviewed Preprint version 1: February 19, 2025
4. Reviewed Preprint version 2: July 3, 2025
5. Version of Record published: August 28, 2025

You can cite all versions using the DOI https://doi.org/10.7554/eLife.104456. This DOI represents all versions, and will always resolve to the latest one.

© 2025, Hironaka and Sumi

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.