Introduction

The search for food in the absence of informative sensory cues is an essential animal behavior¹. A foraging strategy that is observed in nearly all animals is the Area Restricted Search (ARS)². In ARS animals randomly forage for food³, but appetitive sensory cues (like an encounter with food) will cause the animal to restrict their search area by reorienting more frequently, thus increasing the likelihood of food encounters. Conversely, when removed from food, animals will decrease their reorientation rate to increase dispersal^1,4. Caenorhabditis elegans and Drosophila melanogaster larva appear to progressively increase their diffusion constant while foraging off of food by decreasing their rates of reorientation⁴. However in separate studies^5,6, individual worms appear to make an abrupt change from a high to low rate of reorientation. This behavior has been described as a switch from a local to global search strategy, that relies on evidence accumulation to trigger the behavioral switch⁶. In López-Cruz et al.,⁵ the foraging behaviors of individual worms were tracked for 45 minutes after being removed from food. As observed previously⁴, the reorientation rate from this study followed an exponential decay (figure 1a). It was reported that roughly half the worms appeared to make a sudden single switch from local to global search as observed in Calhoun et al⁶, however the other half appeared to produce no discernable change in search strategy, or exhibited multiple switches (figure 1 b-d). Whether or not a worm performed a single decision was defined in Lopez-Cruz et al⁵ by fitting individual reorientation data to two lines using the MATLAB function findchangepts (figure 1b)⁷. This function divides each trace into two regions that are defined by minimizing the sum of the residual squared error of two local linear regressions. The location of the transition point is varied until the total residual error attains a minimum (figure 1b). A large change in slope indicates a sudden change in reorientation rate, and the intersect between the two lines determins the decision-time⁵ (figure 1b). However, the resulting distributions of slope-differences (s1-s2) and transition times for all of the experimental data are continuous, with no clear boundary between deciders and non-deciders (figure 1e).

Figure 1:

Why do some worms appear to make a decision, while others do not? In aggregate, the reorientation rate decays (figure 1a). Despite the population average conforming to a gradual decay, individual trajectories produce a wide diversity of trajectories which sometimes conform to an apparent drop in reorientation rate (figure 1b), while others do not (figure 1c). If the worms are executing a decision, this would seem to indicate only a fraction of the worms decide to switch from local to global foraging strategies, while others perform an alternative strategy.

A central challenge with defining these possible state transitions is the stochastic nature of the behavior itself. Reorientations are random and follow Poisson statistics^8–14. This is very similar to the temporal behavior of individual molecules in solution¹⁵. Individual molecules defined by the same reaction kinetics can stochastically produce long or short time intervals between reaction events, not because one molecule is inherently faster than the other, but because they are stochastically sampling from the same probability distribution. The diversity of times between individual worm reorientations is very similar to this. The exponentially decaying reorientation rate emerges from the average of trajectories that do not necessarily conform to this curve individually. However, even those that produce abrupt switches in reorientation rates could still emerge from a simple exponential decay strategy. Since reorientations occur stochastically, the abrupt changes in reorientation rates could simply be the result of stochastic sampling of an underlying decay phenomenon¹⁶.

Results

We tested this hypothesis by modeling individual worms by stochastic sampling of a decaying reorientation rate with the Gillespie algorithm (figure 2a), a common strategy used to model the kinetics of individual molecules¹⁷. With this strategy, the time between chemical events is modeled by randomly sampling from the time-interval distributions defined by the reaction rates. Although the algorithm was originally developed to model discrete molecular events based on known kinetic parameters, it can be used to generate time trajectories for any discrete events when the kinetics are known. A behavioral example of this is the Lotka-Voltera predator-prey competition model where predator and prey populations fluctuate out of phase due to predation. Stochastic fluctuations of predator and prey populations can be modeled using the Gillespie algorithm^18–20.

Figure 2:

In our model, there are two rates drawn from fitting an exponential decay function to the experimental average reorientation rate (figure 2b):

Where Ω is reorientations α is the initial reorientation rate at t=0, and γ is the reorientation decay rate. To model both the reorientations, and the decay in the reorientation rate, we modeled two kinetic processes, a reorientation rate (a₁) and a decaying rate of parameter M (a₂) that the reorientation rate relies on:

therefore:

where M₀ is the amount of M at t = 0. In this model, the reorientation rate (Equation 2) is a function of a factor M which decays in time (Equation 3), thus producing the decay in the reorientation rate. The initial reorientation rate at t = 0 is α.

To stochastically model both processes, we started with an initial condition of 0 reorientations (Ω₀ = 0) and a large value of M to ensure a smooth decay rate (M₀ = 1000). Since two kinetic events are occurring (the loss of M and the gain of Ω), the algorithm simultaneously generates two time series for M and Ω, by randomly deciding when an event will occur, and the identity of that event. Two random numbers (r₁ and r₂) are drawn from a uniform distribution in [0:1]. One random number (r₁) is used to assign the time of an event (τ) as a combined rate of both rate processes:

where:

A second random number (r₂) is used to calculate which event (j) (1: Ω = Ω + 1 or 2: M = M − 1) occurs based on the relative probabilities of either event:

For example, if the reaction rate a₁ is 10% of a₀ at time t, then there is a 10% chance that the resulting stoichiometry is due to this reaction. In this way, the Gillespie algorithm enables the simultaneous modeling of numerous reactions in parallel with machine-precision time resolution. To ensure the reorientation rate itself is a smooth decay, a large value of M is used to approximate a continuous function. The underlying assumptions are:

1. Reorientations are stochastic^8–14.

2. All the worms experience a common decay of a signaling factor (M) (Equation 3) that influences the reorientation rate (Equation 2).

In our approach, we fit the exponential curve in Eq. 1 to the average reorientation rate from the experimental data from López-Cruz⁵ (figure 2b), and then used these parameters (α and γ) to model 1,631 individual worms (the same number of animals in the experimental data). Starting rates were drawn from a distribution centered at α to match the observed spread of starting rates observed in the experimental data (figure 1a).

In silico reorientation curves generated with this approach produced single-transition and multi/absent transition curves, as observed experimentally (figure 2c). When plotted along with the experimental data, the in silico data produced a distribution of linear regression parameters comparable to the experimental worms (figure 2d). The simulation was able to produce individual trajectories that demonstrated switching behavior, despite the lack of a switching mechanism in the model (Equations 2 and 3). Furthermore, our model demonstrated a continuum of switching to non-switching behavior that was observed in experimental results (figure 2d).

The experimental transition distribution was shifted slightly earlier, however it is important to note that the experimental data were drawn from experiments where 10–15 animals were tracked together. Collisions occurred, but were excluded. However, in addition to touch, animals also reorient in response to pheromones, the presence of which may contribute to experimental observations^21,22. The model also deviated slightly from the experimental data for traces with late inflection points (<3% of experimental data, figure 2d). Again, since the experiments were performed with groups of animals, encounters with pheromone tracks at later time points as the animals dispersed could have potentially altered behavior^21,22.

This exception aside, modeling worm foraging behavior with a simple exponential decay of reorientations was sufficient to capture most experimentally observed dynamics, both switch and non-switch-like. Sudden switches between fast and slow reorientation rates were observed in both the experimental and modeled data (figure 2e), despite the model not relying on a switch strategy, and the decay in M being continuous.

An alternative search strategy is the Lévy walk, where the lengths of time between reorientations (τ) follow a bounded power-law distribution^23,24.

The variables a and b represent the boundaries of shortest and longest observed search lengths, respectively. In this search strategy, bouts of short path lengths randomly occur, independent of sensory information, or internal switching. While this approach also can produce sudden switches in reorientation rates, the average reorientation rate is not time-dependent, unlike what is observed in C. elegans (figure 1a). An advantage of a Lévy walk is that multiple length-scales can be sampled. Characteristically, when the distribution of path lengths is plotted on a log-log scale, the distribution is linear (figure 3) and “scale-free” because there is no characteristic timescale. This is unlike exponential distributions from random walks (cyan through magenta traces, figure 3) where the timescale is defined by the rate constant α.

Figure 3:

There has been considerable debate about the existence of Lévy walks in foraging behavior^25–28, in part because observed search strategies seem to sample multiple timescales, but not to the extent of a true Lévy walk. When these foraging data are compared to only a single exponential or Lévy distribution using maximum likelihood analysis, the results tend to favor the exponential model^27,28, but not exclusively. Observed walk length distributions often have heavier tails than single-exponential distributions, but are not sufficiently sampled to match a Lévy walk²⁸. The distributions observed in this study from the experimental and model animals fall somewhere in between a Lévy walk and single exponential distribution (figure 3). Since the reorientation rate constant decays in time, the resulting distribution has a broader sampling of timescales, similar to a Lévy walk, but bounded by the initial and final observable rates. Since the reorientation rate decays in time, the breadth of this distribution increases with the total search length. This allows the animal to sample from a broader distribution of timescales than simply having one, or even two, reorientation rates.

Discussion

It is important to emphasize that the model was not explicitly designed to match the sudden changes in reorientation rates observed in the experimental data. Kinetic parameters were simply chosen to match the average population behavior. Sudden changes in reorientation rate were not due to sudden changes in the underlying model; stochastic behavior naturally produces sudden bunching of random events. Even if the reorientation rate is set to a constant value, stochastic sampling will still produce sudden changes in rate, even though the rate has no time-dependence (figure 2f).

The lack of a decision simplifies the dispersal strategy for C. elegans. Rather than relying on the accumulation of evidence to make a discrete decision, the worm relies on a decaying signal in the absence of food that drives the reorientation rate. This strategy increases the diffusion constant of the worm, and ensures a more efficient search strategy to find food⁴. The stochastic nature of this search is consistent with prior characterizations of worm behavior and neuronal dynamics^8–14, and implies that any individual worm would exhibit all the variability of the population if allowed to perform this task multiple times. We consider this to be a null model: simple stochastic models should be sufficient to explain observably stochastic behaviors.

The decay in M can be considered the memory timescale of the last food encounter. With an observed γ of 0.07 min^-1, this would mean a memory τ_1/2 of ~10 minutes. In principle M should rise on the presence of food to increase the reorientation rate, on a timescale that is not addressed here. The reliance of the reorientation rate on a decaying factor (M) implies that the rate is influenced by a signaling factor which decays in time. It is well established that the reorientation rate is influenced by numerous signaling factors, and is the basis of the biased random walk that drives taxis in shallow gradients^8–14. The physical basis of M can come from multiple sources. The loss of sensory stimuli alters metabotropic glutamate signaling from sensory neurons which in turn modify the kinetics of the motor network⁵. Altered ionotropic glutamate signaling and dopamine release also influence foraging kinetics²⁹, as well as neuropeptides³⁰. The signaling factor M could be the result of one or all of these factor concentrations decaying in time. Further work will be needed to reveal how the kinetics of reorientations emerges explicitly from underlying signaling kinetics.

Data availability

Experimental data were curated from López-Cruz et al⁵.

Acknowledgements

We thank A. López-Cruz and C. Bargmann, for sharing their experimental data, and also thank them, A. Samuel and members of the Gordus lab for helpful discussions and comments on the manuscript. A.G. acknowledges funding from NIH (R35GM124883).

Additional information

Source Code availability

Analysis was performed using custom written code in MATLAB which can be found in this repository: https://github.com/GordusLab/Margolis_foraging

Author Contributions

A.M. and A.G. conceived the research plan, analyzed the data, and wrote the paper.