A Spiking Neural Network for Autonomous Search and Contour Tracking Inspired by C. elegans Chemotaxis and the Lévy Walk

Abstract

Caenorhabditis elegans exhibits sophisticated chemotaxis behavior through two parallel strategies, klinokinesis and klinotaxis, executed entirely by a small nervous circuit. It is therefore suitable for inspiring fast and energy-efficient solutions for autonomous navigation. As a random search strategy, the Lévy walk is optimal for diverse animals when foraging without external chemical cues. In this study, by combining these biological strategies for the first time, we propose a spiking neural network model for search and contour tracking of specific concentrations of environmental variables. Specifically, we first design a klinotaxis module using spiking neurons. This module works in conjunction with a klinokinesis module, allowing rapid searches for the concentration setpoint and subsequent contour tracking with small deviations. Second, we build a random exploration module. It generates a Lévy walk in the absence of concentration gradients, increasing the chance of encountering gradients. Third, considering local extrema traps, we develop a termination module combined with an escape module to initiate or terminate the escape in a timely manner. Experimental results demonstrate that the proposed model integrating these modules can switch strategies autonomously according to the information from a single sensor and control steering through output spikes, enabling the model worm to efficiently navigate across various scenarios.

1. Introduction

Chemotaxis is a spatial orientation behavior of animals toward chemical cues they deem attractive or away from surroundings they find repellent. Biological chemotaxis strategies are important inspirational resources for autonomous navigation algorithms [1, 2] that control the motion of robots designed for specific tasks, such as locating hazardous chemical leaks, measuring radiation, and environmental monitoring.

Caenorhabditis elegans (C. elegans) exhibits sophisticated chemotaxis toward numerous compounds, including salt (NaCl) [3], with a small-scale and entirely mapped nervous system of only 302 neurons [4, 5], making it an ideal model organism. It is typically attracted to higher salt concentrations for foraging [3], but also exhibits experiential dependence [6], navigating up or down NaCl gradients (positive or negative chemotaxis, respectively) to reach the previously grown salt concentration (setpoint). Furthermore, C. elegans can orient to gradients that vary in intensity by four orders of magnitude [3]. This behavior employs two distinct strategies: klinokinesis [7], a biased random walk, in which reorientation is governed by increasing the frequency of sharp turns away from the desired concentration, and klinotaxis [8], in which the direction of movement is gradually adjusted toward a steeper gradient. Intriguingly, klinotaxis requires the concentration difference between two sides of the body, whereas the salt sensory neurons ASEL and ASER [9] respectively sense the increase and decrease in NaCl at a single point due to proximity [10]. Therefore, chemotaxis in C. elegans is a desirable biological prototype for autonomous navigation models. In particular, klinotaxis enables navigation decisions by comparing left and right concentrations with a single concentration sensor; this makes sense for highly resource-constrained robots [11, 12] that are too small to carry two sensors separated by a certain distance.

Several navigation models based on C. elegans chemotaxis already exist. Morse et al. [13] used a network controller inspired by C. elegans chemotaxis to perform phototaxis in a robot and extracted the underlying computation rule [14]. Dunn et al. [15] simulated klinokinesis behavior using a chemotaxis network model. Deng et al. [16, 17] utilized dynamic neural networks to mimic salt attraction and toxin avoidance and incorporated body undulatory models to produce worm-like robot prototypes. However, the navigation processes of the above models do not actually follow klinotaxis. Izquierdo et al. [18, 19] evolved minimal klinotaxis models and proposed a mechanism underlying positive klinotaxis, whereby steering responses depend on the internal state at the time of the stimulus during sinusoidal locomotion. Chen et al. [20] subsequently identified the same mechanism in a model based on the C. elegans connectome. These studies [18-20] have dual biological and engineering implications, providing the first plausible mechanism for klinotaxis and theoretical support for applying klinotaxis to navigational control.

Meanwhile, recent studies [21, 22] have shown that spiking neural networks (SNNs) [23] are beneficial for navigational control applications. SNNs process information in an event-based manner using spikes, which holds the potential for fast and energy-efficient computing and makes them suitable for hardware implementation [24, 25]. These advantages motivate the transition of navigation models inspired by C. elegans chemotaxis to SNNs. Costalago-Merulo et al. [26] combined an SNN based on the C. elegans connectome with a physical engine to implement klinotaxis, but with two sensory neurons sensing left and right concentrations respectively. Notably, all the aforementioned models achieve positive or negative chemotaxis in isolation, not towards a setpoint. By contrast, Santurkar et al. [27] developed a simple SNN-based klinokinesis circuit to search and track the concentration setpoint. By designing spiking rate-coded logic gates, Shula et al. [28] improved the circuit in [27] and further developed an escape circuit using an XOR logic gate for the local extrema problem; nevertheless, the escape duration is fixed by injecting an external inhibitory current for a specific duration. Kishore et al. [29] applied the klinokinesis circuit [27] to Loihi [24], a neuromorphic chip. Shukla et al. [30] extended the state-dependent mechanism of positive klinotaxis [18, 19] to negative klinotaxis and incorporated a klinotaxis module into the aforementioned model [28]; nevertheless, the module processes the received spikes mathematically rather than with an SNN circuit.

Chemotaxis is effective in cases where gradients can be detected but fails in the absence of a gradient. In practice, robots sometimes must “blindly” explore the environment without an initial gradient. The Lévy walk, as a specialized random search strategy, is optimal for searching sparse resources without external cues and has been observed in humans [31] and other animals [32-34]. It comprises a cluster of short steps intermittently interspersed with longer movements, allowing for farther travel from the initial position [35] and thus increasing the encounter rate with targets. Therefore, combining C. elegans chemotaxis with the Lévy walk is expected to achieve efficient searches in complex scenarios with and without gradients.

Based on the above premises, we propose an SNN model jointly inspired by chemotaxis in C. elegans and the Lévy walk. It aims to search and track concentration setpoints using a single sensor across complex scenarios with and without gradients and with local extrema. The main contributions are as follows:

(1) A klinotaxis module using spiking neurons is designed based on the existing non-SNN klinotaxis model. This module controls steering with spikes and allows rapid searches along almost the steepest gradients.

(2) A random exploration module is built, which generates Lévy-walk motion when no gradient is detected to increase the encounter rate with gradients. This is the first realization of Lévy walks using an SNN model.

(3) The escape circuit is improved, and a termination module is further constructed. This two module can make decisions about the escape and termination of the escape in time to leave the trap and switch to other appropriate strategies.

(4) All modules are integrated seamlessly into the integrative model through synaptic connections. Experimental results demonstrate that the model can autonomously switch navigation modes according to sensory information; thus, it can control the model worm to quickly encounter gradients, escape from traps, reach the setpoint with a short path, and eventually track the setpoint contour with a small deviation.

2. Methods

2.1 Overall Architecture and Neuron Model

The overall model architecture, as shown in Fig. 1, consists of an information preprocessing module and an SNN. The information preprocessing module generates current inputs to the SNN according to the concentration information collected by a single sensor and the concentration setpoint. The SNN contains six sub-network modules with different functions, among which the sensory module encodes sensory-related inputs and the other modules control various navigation strategies with spike trains.

Fig. 1. Overall model architecture. The gray part is the information preprocessing module, where squares represent functional units. The other colored parts are SNN modules, where circles represent neurons. The six color-filled circles are output neurons that control steering. The black, blue, and red lines represent the external injection currents, excitatory and inhibitory synaptic connections of the neurons, respectively.

All neurons in SNN are modeled as leaky integrate-and-fire (LIF) neurons with current-based synapses [36]. The dynamic equation of neuron i is given as follows.

\(\begin{aligned}\tau_{m, i} \frac{d V_{i}(t)}{d t}=-V_{i}(t)+I_{i}(t)\end{aligned}\),       (1)

where τ_m is the membrane time constant. Here V( t ) denotes the membrane potential and I( t ) is the input current, as shown below:

\(\begin{aligned}I_{i}(t)=\sum_{j} \sum_{t_{j}^{k}<t} w_{i j} \cdot \exp \left(-\frac{t-t_{j}^{k}}{\tau_{s, i j}}\right)+I_{i n j, i}(t) \end{aligned}\).       (2)

The first term represents synaptic currents from other neurons, where τ_s is the synaptic time constant and w_ij is the synaptic strength from neuron j to neuron i. Here w_ij > 0 (w_ij < 0) corresponds to the excitatory (inhibitory) connection and the t^k_j is the time when neuron j emits the k- th spike. I_inj( t ) represents the external injection current.

When the input current of neuron i pushes its membrane potential above the threshold potential V_th,i, the neuron activates, fires a spike, and clamps its membrane potential to the reset value (here taken to be zero). After spiking, the neuron becomes refractory for a time period τ_ref,i; that is, V_i(t)= 0 for t^k_i < t ≤ t^k_i +τ_ref,i (k =1, 2, …). Then, its output spike train is as follows.

\(\begin{aligned}S_{i}(t)=\left\{\begin{array}{ll}1 & t \in \Gamma_{i} \\ 0 & \text { otherwise }\end{array},\right.\end{aligned}\)       (3)

where Γ_i ={t¹_i, t²_i, t³_i, …} is the set of firing times of neuron i.

The spike trains of six output neurons from three modules control steering. We assume a model worm whose instantaneous angular direction at time t is θ(t).

\(\begin{aligned}\frac{d \theta(t)}{d t}=\sum_{i} \omega_{i} \cdot\left(S_{i r}(t)-S_{i l}(t)\right).\end{aligned}\)       (4)

The subscript i ∈ {p, w, r} corresponds to the klinokinesis, klinotaxis, or random exploration modules, respectively. Here ω > 0 is the turning bias caused by a spike of the corresponding output neuron. Studies [8, 14] have shown that C. elegans moves forward at an approximately constant speed. For simplicity, we make the model worm move at a constant speed v, and its position (x(t), y(t)) at time t is updated as follows:

\(\begin{aligned}\left(\frac{d x(t)}{d t}, \frac{d y(t)}{d t}\right)=(v \cdot \cos \theta(t), v \cdot \sin \theta(t)).\end{aligned}\)       (5)

Fig. 2 shows the flow diagram for the navigation algorithm implemented by our model, which indicates that the model needs to autonomously switch between different navigation modes according to the currently sensed scene. Subsequently, we depict each module in detail to show how the modules work collaboratively to implement this navigation process; the descriptions of the information preprocessing module will be spread across the sections involving it.

Fig. 2. Flow diagram for the navigation algorithm implemented by the SNN model

2.1 Sensory module

The sensory module contains two setpoint-detection neurons (N_s1 and N_s2) and two gradient-detection neurons (N_s3 and N_s4), as shown in the pink part of Fig. 1. N_s1 and N_s2 receive currents I_s1(t) and I_s2(t), respectively, which are generated by the M1 unit of the information preprocessing module by comparing the instantaneous concentration C(t) from the sensor with the setpoint C_T.

I_s1(t) = I_s0ㆍH(C(t) - C_T + ε₁),       (6)

I_s2(t) = I_s0ㆍH(C_T + ε₁- C(t)),       (7)

where I_s0 is a positive constant and ε₁ > 0 is a small preemptive degree that can alleviate the latency of the steering responses when the setpoint is reached. H(ㆍ) is the Heaviside function. N_s1/N_s2 emits spikes at a fixed frequency when C(t) is below/above C_T, as shown in Fig. 4a.

N_s3 and N_s4 receive currents I_s3(t) and I_s4(t) generated by the M2 unit, respectively. The M2 unit first estimates the temporal gradient via a derivative operator applied to the recent historical concentrations, i.e., C(t) - C_avg(t), where C_avg(t) is the average concentration over the past second. In a previous model [30], an adaptive dynamic scaling based on historical concentrations was adopted to enable gradient-detection neurons to operate over a wide range of concentrations. However, considering two scenarios where the gradients are similar but concentrations are greatly different, the input currents are different and vary widely. Like the sensory neurons ASEL and ASER of real C. elegans, the inputs to gradient-detection neurons should depend on gradients rather than absolute concentrations. We therefore introduce an adaptive dynamic scaling based on historical gradients, which makes the input currents dependent of gradients and independent of absolute concentrations and allows neurons to operate effectively over a wide range of gradients. The current is as follows.

\(\begin{aligned}I_{e s}(t)=\frac{a_{1} \cdot\left[C(t)-C_{a v g}(t)\right]}{1+a_{2} \cdot d C_{a v g}(t)},\end{aligned}\)       (8)

I_s3(t) = max(I_es(t), 0),       (9)

I_s4(t) = max(-I_es(t), 0),       (10)

where a₁ and a₂ are positive constants, and dC_avg(t) is the average |C(t) - C_avg(t)| over the past period. Therefore, N_s3/N_s4 receives the current and emits spikes when the gradient is positive/negative; moreover, the spiking frequency is proportional to the gradient magnitude and the maximum frequency is limited by the refractory period. Fig. 4a shows the spiking responses in a concentration profile.

2.3 Klinotaxis Module

According to the state-dependent mechanism of positive klinotaxis [18, 19] and that of negative klinotaxis extended by [30], the klinotaxis strategy relies on the alternating saturation of outputs that control left and right turns at different phases during sinusoidal locomotion. Due to the threshold potential and refractory period, the spiking frequency of spiking neurons naturally exhibits a saturated nonlinear relationship with the input current, with a sensitive region and two saturated regions, as shown in Fig. 3a. Exploiting this characteristic, we design a klinotaxis module with two interneurons (N_w1 and N_w2) and two output neurons (N_wr and N_wl), as shown in the blue part of Fig. 1.

Fig. 3. Neuronal dynamics of the klinotaxis module and resulting steering responses. (a) Spiking frequency vs. input current of interneurons in the absence of sensory inputs. (b)-(e) Changes in the responses of output neurons (bottom) and motion trajectories (top) emitted by the concentration increase (red) or decrease (green) when the bias current is b₁ (i.e., positive search scenario) (b, c) or b₂ (i.e., negative search scenario) (d, e).

To generate sinusoidal motion, N_w1 and N_w2 respectively receive oscillating currents with opposite phases generated by the M4 unit.

\(\begin{aligned}I_{o s c 1}(t)=-I_{o s c 2}(t)=I_{o s c 0} \cdot \sin \left(\frac{2 \pi}{T_{o s c}} \cdot t\right),\end{aligned}\)       (11)

where I_osc0 > 0 is the maximum oscillation amplitude of the current, and T_osc = 4.2 s is the oscillation period, which is consistent with that of real C. elegans. Meanwhile, N_w1 and N_w2 receive the same bias current I_b(t) generated by the M3 unit according to the relationship between C(t) and C_T.

\(\begin{aligned}I_{b}(t)=\left\{\begin{array}{ll}b_{1} & C_{T}-C(t)>\varepsilon_{2} \\ b_{3} & \left|C(t)-C_{T}\right|<\varepsilon_{2} \\ b_{2} & C(t)-C_{T}>\varepsilon_{2}\end{array}\right.,\end{aligned}\)       (12)

where ε₂ > 0 is a small constant, and b₁ > b₂ > b₃. As shown in Fig. 3a, the oscillating input currents of N_w1 and N_w2 are shifted to different regions by superposing I_b(t) ; the case for b₃ is not shown because b₃ shifts currents to the left into the fully saturated region regardless of the oscillating phase. By doing so, the model will autonomously switch between positive klinotaxis for C(t) < C_T - ε₂, negative klinotaxis for C(t) > C_T + ε₂, and the close of oscillation and klinotaxis for |C(t)- C_T| < ε₂. In addition, N_w1 and N_w2 receive the same sensory inputs, with excitatory and inhibitory connections from N_s3 and N_s4 (Fig. 1); that is, w_w1,s3 = w_w2,s3 = -w_w1,s4 = -w_w2,s4 > 0_.

The output neurons N_wr and N_wl receive the difference between the outputs of N_w1 and N_w2 with antagonistic synaptic effects; that is, w_wr,w1 = -w_wl,w2 = -w_wl,w1 = w_wl,w2 > 0. When the spiking frequency of N_w1 is greater than that of N_w2, N_wr receives the positive input current and emits spikes; the reverse is true for N_wl. The spike trains of N_wr and N_wl respectively cause right and left turns.

The dynamics of N_w1 and N_w2 in the absence of sensory input are shown in Fig. 3a. Influenced by the anti-phase sinusoidal currents, the spiking frequency of N_w1 is greater than that of N_w2 during the first half of the oscillatory cycle, thereby N_wr being activated; N_wl is activated during the other half of the cycle. Consequently, the model worm moves forward in a sinusoidal wave pattern when it is not near the setpoint (I_b(t) ≠ b₃), as shown in the blue curves in Fig. 3b-e.

Next, we describe how this module generates steering responses that produce klinotaxis in the presence of sensory inputs. For concentrations less than C_T (i.e., positive search scenario), spiking frequencies of N_w1 and N_w2 alternately enter the upper saturation region (green and blue curves in Fig. 3a) due to the large positive bias current b₁. Fig. 3b and 3c show the changes in the spiking responses of N_wr and N_wl and the steering responses caused by gradients at different phases. When a positive/negative gradient is sensed at a certain phase, N_s3/N_s4 emits spikes, thus increasing/decreasing the input currents of N_w1 and N_w2 over a subsequent period. The spiking frequency of the neuron entering the saturated region remains unchanged or changes slightly, whereas the frequency of the other neuron in the sensitive region increases/decreases substantially. That is, the positive/negative gradient causes the subsequent output difference between N_w1 and N_w2 to decrease/increase. Accordingly, the number of spikes of N_wr or N_wl decreases/increases, and thus the subsequent turning angle of the oscillating motion decreases/increases. As such, the translational direction is adjusted to the same/opposite side of the instantaneous direction of that phase. This is effective because the presence of a positive/negative gradient at a phase point implies that the instantaneous direction at that point is the direction of increasing/decreasing concentration. Consequently, the klinotaxis module gradually adjusts the model worm towards a higher concentration (i.e., positive klinotaxis).

For concentrations greater than C_T (i.e., the negative search scenario), N_w1 and N_w2 alternately enter the lower saturation region (the purple and red curves in Fig. 3a) due to the near-zero bias current b₂. Fig. 3d and 3e show the responses in this case. Spikes in N_s3/N_s4 (positive/negative gradient) result in an almost constant output of the neuron in the saturated region, along with an increased/decreased output of the neuron in the sensitive region. The subsequent output difference between N_w1 and N_w2 increases/decreases. Accordingly, the number of spikes of N_wr or N_wl increases/decreases, resulting in subsequent increases/decreases in the turning angle. As such, the translational direction is adjusted to opposite/same side of the instantaneous direction of that phase, that is, the direction of the decreasing concentration. Consequently, this module gradually directs the model worm towards a lower concentration (i.e., negative klinotaxis).

For concentrations close to C_T, the large negative bias current b₃ inhibits N_w1 and N_w2, causing them to fully enter the saturation region. They are not activated regardless of the oscillating current and sensory input. Therefore, N_wr and N_wl always lack an input current and do not emit spikes. In this case, the oscillating motion and klinotaxis are both eliminated, and only the klinokinesis strategy is used to track the contour of the setpoint.

Additionally, we set up inhibitory connections from the escape module to N_wr and N_wl; thus, once the model worm gets stuck in a local trap and begins to escape, the klinotaxis module is turned off.

It can be observed from Fig. 3a that the input current range to which interneurons N_w1 and N_w2 are sensitive is approximately one. This sensitive region should be fully utilized to simultaneously ensure sufficient motion oscillations and the sensitivity to the sensory input; therefore, in our simulation, the maximum afferent current from gradient-detection neurons is 0.4 by setting synaptic strengths, and I_osc0 = 0.6. Additionally, because the input-frequency curve (the black curve in Fig. 3a) is steeper near the lower saturation region, the spike frequencies of N_w1 and N_w2 during negative klinotaxis (I_b(t) = b₂) vary more widely with oscillating currents than that during positive klinotaxis (I_b(t) = b₁), resulting in more intense oscillating motion (compare Fig. 3b and Fig. 3d). The amplitude of the frequency variation with the oscillating current can be reduced by decreasing b₂; however, if b₂ is too small, the oscillating motion cannot be ensured for all sensory inputs (i.e., when the absolute value of the negative gradient is large, N_w1/N_w2 may be always in the lower saturation region). Therefore, we set b₂ slightly less than the sensitive threshold current.

2.4 Klinokinesis Module

The klinokinesis module follows the design in [28], additionally adding inhibitory connections from the escape module. This module contains only two motor neurons (N_pr and N_pl), as shown in the green part of Fig. 1. These two neurons and sensory neurons constitute two “AND” analog logic gates in order to realize klinokinesis. Let O_i denote the spiking frequency of neuron i. Then,

O_pr = (￢O_e5)∧(O_s1 ∧ O_s3),       (13)

O_pl = (￢O_e5)∧(O_s2 ∧ O_s4),       (14)

where ¬ and ∧ denote logical operators “NOT” and “AND”, respectively. Equation (13) indicates that N_pr emits spikes when N_s1 and N_s3 simultaneously emit spikes and N_e5 does not; its spiking frequency is proportional to that of N_s3. A similar understanding holds for (14). That is, N_pr or N_pl is only activated when the model worm is not in the escape mode and is moving away from the setpoint. Specifically, if C(t) > C_T and the gradient is positive, N_pr emits spikes, leading to right turns; if C(t) < C_T and the gradient is negative, N_pl emits spikes, leading to left turns. Moreover, the greater the deviation between the motion direction and the setpoint direction, the higher the spiking frequency of N_pr or N_pl, resulting in a stronger directional correction. Fig. 4a shows the spiking responses for implementing klinokinesis in an example concentration profile.

By setting threshold potentials and afferent synaptic strengths of N_pr and N_pl, the fixed spiking frequency of N_s1/N_s2 makes the voltage of N_pr/N_pl close to the threshold potential, but not enough to activate N_pr/N_pl; if at the same time N_s3/N_s4 emits spikes, N_pr/N_pl is activated and its spiking frequency is proportional to that of N_s3/N_s4. Furthermore, the inhibitory current from N_e5 must be large enough to inhibit N_pr/N_pl regardless of sensory inputs; once the escape mode is enabled, N_pr and N_pl are inhibited and the klinokinesis module is automatically disabled.

Fig. 4b shows the spiking responses when the model worm tracks the setpoint contour. Once the setpoint is reached and overshot, either N_pr or N_pl is activated, causing turns towards the setpoint; subsequently, N_pr and N_pl emit spikes alternately, allowing oscillating motion along the setpoint contour. Additionally, it can be observed from Fig. 4 that ε₁ in (6) and (7) causes N_s1 and N_s2 to be activated before reaching the setpoint; this reduces the response latency of N_pr and N_pl around the setpoint, thereby reducing the deviation during contour tracking.

Fig. 4. Spiking responses of sensory and klinokinesis modules. (a) Responses in an example concentration profile. The profile is divided into four regions according to the relationship between C(t) and C_T and the algebraic sign of gradients. (b) Responses during contour tracking.

2.5 Random Exploration Module

The Lévy walk is a specialized random walk, comprising mostly short move steps (distance traveled before pausing or changing direction) interspersed with rarer longer steps. The move-step length is distributed according to a power-law function with a heavy tail: P(l)~l^−µ (~ means “distributed as”) with 1<µ≤,3 where l is the move-step length and µ is the Lévy exponent. A random exploration module is designed to generate Lévy-walk motion in a flat-concentration scenario, as shown in the yellow part of Fig. 1. The module consists of three interneurons (N_r1, N_r2, and N_r3) and two motor neurons (N_rr and N_rl).

N_r1 and N_r2 respectively receive currents I_a1(t) and I_a2(t) from the M5 unit. I_a1(t) = max(I_a(t), 0) and I_a2(t) = max(-I_a(t), 0), where I_a(t) is a current with an update interval of one second. At each update, I_a(t)=0 follows a probability distribution with a power-law tail:

p(I_a(t > t_n) = 0) = (t - t_n)^-μ,       (15)

where t_n is the last time I_a(t) became non-zero (i.e., t_n is updated once every time I_a(t) changes from zero to non-zero), and µ = 2 which has been shown to be the theoretically optimal exponent for the highest searching efficiency [32, 33]. As shown in Fig. 5a, the non-zero I_a(t) will be updated to zero with a probability of one in the next second; the longer I_a(t) is zero, the higher the probability of becoming non-zero. The positive/negative I_a(t) activates N_r1/N_r2, emitting spikes whose frequency is proportional to the magnitude of I_a(t).

Fig. 5. Response dynamics of the random exploration module. (a) The probability that I_a(t) is zero per update. When t = 1 s and t = 7 s, I_a(t) becomes non-zero (t_n changes from 1 to 7). (b) Responses dynamics in an example concentration profile.

N_r3 has excitatory connections from N_wr and N_wl and inhibitory connections from N_s3 and N_s4, which determines whether to turn on the random exploration mode.

O_r3 = ￢(O_s3 ∨ O_s4 )∧ (O_wr ∨ O_wl),       (16)

where ∨ denotes the logical operator “OR.” Equation (16) indicates that N_r3 is activated if there is no gradient (i.e., neither N_s3 nor N_s4 fires) and the setpoint has not been reached (i.e., N_wr or N_wl fires). To form this function, N_r3 has a small threshold potential and large afferent synaptic strengths from N_wr and N_wl; thus, a single spike from N_wr or N_wl can activate N_r3 and cause it to emit spikes at the maximum frequency. Furthermore, inhibitory connections from N_s3 and N_s4 are large enough to inhibit N_r3 in the presence of any concentration gradient.

The output neurons perform analog “AND” operations on spike trains of interneurons.

O_{rr =}  O_r1 ∧ O_r3 ,       (17)

O_{r1 =} O_r2 ∧ O_r3 .       (18)

The spiking frequencies of N_rr and N_rl is proportional to the magnitude of I_a(t) through the transmission of N_r1 and N_r2. Their spike trains control right and left turns, respectively.

The output spikes caused by non-zero I_a(t) interrupt the move step and reorient the model worm. The study [33] has found that the turning angle distribution for Drosophila larva paths fitted by truncated power-laws is a near uniform distribution. Therefore, in our simulation, the updated non-zero I_a(t) are random values uniformly distributed within a certain range, so that the number of output spikes also varies nearly uniformly. Consequently, in the absence of gradients, this module autonomously activates, inducing the model worm to perform Lévy-walk motion to randomly explore the environment; Lévy-walk motion involves a power-law distribution of the move-step length and a near uniform distribution of the turning angle. Fig. 5b shows the response dynamics of this module in an example concentration profile.

2.6 Escape and Termination Modules

When stuck in a local maximum/minimum during positive/negative klinotaxis, concentrations are below/above the setpoint. Once the model worm moves away from the “peak”/“valley”, the detected concentration decreases/increases, and N_pl/N_pr emits spikes making it turn back; repeatedly, the model worm will linger at the “peak”/“valley”. In this case, N_pl/N_p fires exclusively in the klinokinesis module. [28] has developed an XOR gate-based escape circuit to determine whether the simulation robot is in a local trap. We adopt this circuit and improve it so that the circuit can interact with other modules without external inhibitory currents, as shown in the red part of Fig. 1.

N_e1 and N_e2 can be activated by a single spike from N_pr and N_pl, respectively, and have long refractory periods. N_e3/N_e4 acts as an integrator, which is activated when N_e1/N_e2 emits three consecutive spikes and N_e2/N_e1 does not. N_e5 performs an “OR” operation on N_e3 and N_e4; a single spike can activate N_e5 and cause it to emit spikes at the maximum frequency. Consequently,

O_e5 = O_pl ⊕ O_pr ,       (19)

where ⊕ denotes the logical operator “XOR.” This indicates that N_e5 is activated only if either N_pr or N_pl emits spikes continuously for a while and the other does not.

There are several improvements compared to the previous escape circuit [28]. First, for N_e5, we added the self-connection, inhibitory connections to other modules, and inhibitory connections from the termination module. In so doing, we ensure that once N_e5 is activated, it continues to emit spikes, inhibiting other navigation modules, so that the model worm moves straight away from the trap until N_e5 receives inhibitory signals from the termination module. In addition, we add a new neuron N_e6, which is activated after N_e5 fires continuously for a certain period due to its large threshold potential and afferent synaptic time constant. N_e6, rather than N_e5, outputs to the termination module to ensure that termination judgments are immune to the oscillating motion prior to the escape mode. Fig. 6a shows the spiking responses of the escape module for a local maximum trap.

Fig. 6. Spiking responses of the escape module for a local maximum trap (a) and spiking responses of termination modules for three different scenarios (b-d). The top diagrams represent concentration profiles. The bottoms are spike patterns of neurons near the moment of escape or termination of escape.

We noticed that when escaping from local maxima/minima, the model worm moves straight along negative/positive gradients. Once it leaves the trap, either the algebraic sign of the gradient changes (Fig. 6b, c) or the gradient disappears (Fig. 6d). Based on this, we develop a termination module that can autonomously decide whether to terminate the escape. The termination module is divided into three portions, as shown in the purple part of Fig. 1.

The first portion, consisting of N_t1 to N_t4, constitutes three “AND” logic gates. It functions in scenarios where gradients exist outside the local maximum trap.

O_t3 = O_t1 ∧ O_s3,       (20)

O_t2 = O_e6 ∧ O_s2,       (21)

O_t4 = O_t2 ∧ O_t3.       (22)

N_t1 is activated by spikes from N_s4 and continues to emit spikes for some time after N_s4 stops spiking due to its large afferent synaptic time constant. Thus, (20) indicates that N_t3 is activated if the gradient changes from negative to positive (condition 1). Equation (21) indicates that N_t2 is activated if the escape is initiated in a local maximum trap (condition 2). Satisfying both conditions means that the model worm has escaped from the trap; then, N_t4 emits spikes to inhibit N_e5 and N_e6, thereby terminating the escape. Fig. 6b shows the spiking responses when this portion is in operation.

To form this function, a single spike from N_s4 can activate N_t1 and cause it to emit spikes at the maximum (fixed) frequency. The fixed spiking frequency of N_t1 makes the voltage of N_t3 close to its threshold potential, but not enough to activate N_t3; if at the same time N_s3 emits a spike, N_t3 will be activated. N_e6 and N_s2 emit spikes at a fixed frequency, resulting in the fixed frequency of N_t2. Similarly, the fixed frequency of N_t2 is not enough to activate N_t4; if at the same time N_t3 emits a spike, N_t4 will be activated.

The second portion, consisting of N_t5 to N_t8, works similarly to the first portion, except that it functions in scenarios where there are gradients outside the local minimum trap. N_t5 is activated by spikes from N_s3 and continues to emit spikes for some time after N_s3 stops spiking.

O_t7 = O_t5 ∧ O_s4,       (23)

O_t6 = O_e6 ∧ O_s1,       (24)

O_t8 = O_t6 ∧ O_t7.       (25)

Fig. 6c shows the spiking responses of this portion.

The third portion, N_t9, operates in scenarios where there is no gradient outside the local trap.

O_t9 = O_e6 ∧ (￢(O_s3 ∨O_s4 ∨ O_t1 ∨O_t5))_.       (26)

If the escape mode is on and no gradient is sensed for a while, N_t9 emits spikes to terminate the escape. Fig. 6d shows the spiking responses for such a scenario.

In conclusion, the termination module can autonomously determine whether the model worm has escaped from traps in various local extremum scenarios and terminate the escape.

3. Experiment Results

We tested the performance of our proposed model in various simulation scenarios, which were 100 mm× 100mm square areas with distinct concentration distributions. The constant speed of the model worm was set to 0.1 mm/s. The ordinary differential equations were solved using Euler integration with a time step of 0.01 s.

3.1 Search and Contour tracking

We measured the search and tracking performance in 20 different scenarios with Gaussian-shaped concentration gradients to test the effectiveness and generalization of our model. The scenarios include positive and negative search scenarios (e.g., Fig. 7c, d), each with different contour shapes and gradient magnitudes (e.g., Fig. 7c, e). The elliptic contour is superimposed by two Gaussian distributions. In each scenario, the model worm started moving with 10 different initial orientations. Furthermore, to verify the effectiveness of the klinotaxis module, we performed the same experiment for the model worm which ignores the output control of the klinotaxis module (i.e., the klinokinesis-only model worm). The search performance was measured by the ratio of the search time spent by the model worm to the search time corresponding to the shortest path (SSR). The tracking performance was measured by the average deviation rate (ADR) from the setpoint during tracking, referring to the metric used in [30].

Fig. 7. Several trajectories of model worm and klinokinesis-only model worm in different scenarios. (a), (b) The klinokinesis-only model worm in positive (a) and negative (b) search scenarios. (c)-(f) The model worm in positive (c, e) and negative (d, f) search scenarios. The black solid lines are trajectories, and the black dotted lines are C_T contour lines. C_T = 55 mM for (a), (c), and (e). C_T = 25 mM for (b), (d), and (f).

\(\begin{aligned}A D R=\frac{1}{T-T_{0}} \int_{T_{0}}^{T} \frac{\left|C(t)-C_{T}\right|}{C_{\max }-C_{\min }} d t\end{aligned}\),       (27)

where T₀ is the time at which the model worm first arrives at C_T , and T is the total movement time. C_max and C_min are the maximum and minimum concentrations in the scenario, respectively.

As can be observed from the trajectories (Fig. 7c, d), regardless of the positive or negative search, the model worm realized klinotaxis and klinokinesis strategies simultaneously and moved forward in a sinusoidal wave pattern. Klinokinesis allowed for quick turns when the model worm was moving in an unfavorable direction and klinotaxis adjusted its direction to a direction with nearly the steepest gradient, ensuring short search trajectories. After the setpoint was reached, klinotaxis and oscillating motion were automatically eliminated and only the klinokinesis strategy was used to alternate between slight left and right turns around the setpoint, thereby tracking the contour. In contrast, the search trajectories of the klinokinesis-only model worm, as shown in Fig. 7a and 7b, were longer and non-oscillating; it turned only when away from the setpoint, and once a favorable gradient was sensed, it moved straight without further adjustment. Table 1 shows the average performance of the two models with different initial orientations in the scenario shown in Fig. 7a or 7b. Compared to the klinokinesis-only model worm, the worm model had significantly shorter search times and comparable ADRs, and performance was independent of the initial orientation (small standard deviation in the SSR). The search time for negative searches was slightly longer than that for positive searches. This was because the spiking frequency of N_w1 and N_w2 varied more with the oscillating current during the negative search (Fig. 3a), resulting in a greater amplitude of the motion oscillation (compare Fig. 7c and 7d).

Table 1. Search and tracking performance of the model worm and the klinokinesis-only model worm in the scenario shown in Fig. 7a or 7b (±s.d.)

Table 2 shows the average performance across all scenarios. The ratio of the average SSR of the model worm to that of the klinokinesis-only model worm was 0.26, that is, the time to reach ratio (TRR) of the two models, which is a metric used in [30]. For the previous non-SNN klinotaxis model [30], the average TRR and ADR were 0.39 and 0.44%, respectively. This demonstrates that our model not only realizes klinotaxis using spiking neurons but also significantly outperforms the previous model in search and tracking performance.

Table 2. Search and tracking performance of the model worm and the klinokinesis-only model worm in all scenarios (±s.d.)

Additionally, we compared the search and tracking performance of the model worm with that of a model worm without adaptive dynamic scaling. The experiments were performed in scenarios with conical gradients, where the gradient is independent of the spatial position and therefore suitable for observing the effect of gradients on performance. The concentration is linearly proportional to the Euclidean distance r from the concentration peak: C(r)∝ gㆍr , where g is the conical gradient coefficient (i.e., the spatial gradient in the steepest gradient direction). As can be seen from Fig. 8, when the concentration changed gently (i.e., g was very small), the search and tracking performance of the non-adaptive model worm decreased sharply; in particular, it cannot search and track when g was less than 0.2. This was because in this case, the input current of the gradient-detection neuron N_s3 or N_s4 was too small for the neuron to be activated, causing severe damage to both klinokinesis and klinotaxis. When the gradient was steep (i.e., g was very large), the search performance of the model worm without dynamic scaling deteriorated. This was because the input current was so large that the spiking frequency of N_s3 or N_s4 was close to saturation, thus the sensitivity of N_s3 or N_s4 to the gradient amplitude being small, affecting the realization of klinotaxis. In this regard, the adaptive dynamic scaling adjusted the sensitivity to gradients according to the recent historical gradients, ensuring the relatively invariant scale of the input currents of the gradient-detection neurons; therefore, the performance of the model worm remained stable over a wide range of gradients.

Fig. 8. Performance comparison of the model worm and the non-adaptive model worm. (a) ADR vs. g . (b) SSR vs. |g| . Each SSR was the average SSR of both positive (g > 0) and negative (g < 0) searches.

3.2 Random Exploration

We placed the model worm in an initial position with a flat concentration (i.e., no gradient) to test the random exploration module. From several trajectories shown in Fig. 9a, groups of short move steps were interspersed with longer move steps between them, analogous to the motion pattern of L é vy walks. Once a gradient was encountered, the model worm automatically switched from L é vy-walk motion to searching along the gradient using klinotaxis and klinokinesis strategies.

Additionally, we compared the exploration efficiency of Lévy-walk motion exhibited in our model with the simple random walk executed in the previous model [27]. For the random-walk model worm, it reoriented once per second at a random turning rate (º/s) uniformly distributed between [-φ_max, φ_max] , and φ_max was consistent with the maximum turning rate of Lévy-walk motion. For the scenario without gradients outside the circle shown in Fig. 9b, the model worm started moving at the same initial position for a maximum of 1,500 s; the experiment was performed 100 times for each strategy. The trajectory of our model worm (Fig. 9b) was more diffuse, while that of the random-walk model worm (Fig. 9c) was more concentrated near the initial position. Table 3 shows the number (m) of times the two types of model worms reached the circle within 1,500 s, the average time (t_random) they taken to reach the circle, and the average shortest distance (D_shortest) from the circle if they failed. By comparison, our model worm had a significantly higher arrival rate in the specified time, a shorter average time, and a smaller average shortest distance. These indicate that the Lévy walk expands the exploration coverage and increases the chance of encountering gradients.

Fig. 9. Random exploration in the absence of the gradient. (a) Several trajectories of the model worm. The white dotted lines and black dotted lines represent the boundaries with gradients and contours of the setpoint, respectively ( C_T = 42 mM ). (b), (c) A trajectory of the Lévy walk (b) or the random walk (c). Inside the circle is the region with gradients. (d) Move-step series ( n = 21,640 ). (e) The move-step time frequency distribution of all move steps. Inset, the log-log plot of move-step frequency versus move-step time. The black solid line is the theoretical distribution of the ideal Lévy walk with µ = 2.

Table 3. Performance of multiple random explorations in scenario in Fig. 9b (n = 100, ± s.d.)

Further, to verify that the random exploration motion indeed conformed to the Lévy walk, a statistical analysis was conducted on all move steps ( n = 21,640 ) of 100 experiments. The move-step time was proportional to the move-step length due to the constant speed. Fig. 9d shows the move-step time series, with an intermittent structure of longer steps. The move-step time frequency distribution of all steps, as shown in Fig. 9e, closely resembled an inverse-square power law with a heavy tail of increasingly longer steps intermittently distributed within the time series, which is typical of ideal Lévy walks. The log-log relationship of the frequency versus move-step time approximated the theoretical distribution with L é vy exponent µ = 2.

3.3 Escape from Local Extremum Traps

To test escape and termination modules, the model worm was placed in various local traps and moved freely. When located in a local maximum trap (Fig. 10a), the model worm first searched the setpoint through oscillating motion up gradients. After hovering around the extreme point for short periods, it entered escape mode, moving straight without oscillation to leave the trap. Escape was terminated at positions where the algebraic sign of gradients changed or gradients disappeared. The model worm then continued to search the setpoint along the gradients or explore randomly to search for the gradient. For local minimum traps, the results were similar except for moving down gradients (Fig. 10b). These indicate that the model can initiate and terminate the escape in a timely manner according to the sensed concentration information, and switch to an appropriate strategy to autonomously complete the navigation task of searching and tracking setpoints.

Fig. 10. Several trajectories in scenarios with local extremum traps. (a) Local maximum and positive search scenario. (b) Local minimum and negative search scenario. The black dotted lines are setpoint contours (C_T = 42 mM (a) or 24 mM (b)). Middle and right, enlarged views of two trajectories on the left, showing terminations of the escape caused by the change in the algebraic sign of the gradient and the gradient disappearance, respectively. The pentacles are initial positions.

4. Conclusion

In this study, we propose a bioinspired and efficient SNN model for search and contour tracking of specific concentrations, which for the first time combines chemotaxis in C. elegans and the Lévy walk. On the one hand, the model fully draws on the advantages of various biological strategies to perform fast searches in scenarios with and without gradients and simultaneously solves the local extrema problem, which may be encountered in processes. A dynamic scaling based on historical gradients is introduced, enabling the model to operate over a wide range of gradients. On the other hand, the model is based on an SNN with a relatively simple architecture and requires only a single sensor. This has the potential to reduce computational costs and power consumption and is beneficial to hardware implementation. We demonstrate the effectiveness and flexibility of the proposed model and its superior performance over previous models across various scenarios. This model can serve as a prototype for robot controllers. We intend to apply the model to an actual robot and validate it in real-world environments in future works.