A neural coding scheme reproducing foraging trajectories

The movement of many animals may follow Lévy patterns. The underlying generating neuronal dynamics of such a behavior is unknown. In this paper we show that a novel discovery of multifractality in winnerless competition (WLC) systems reveals a potential encoding mechanism that is translatable into two dimensional superdiffusive Lévy movements. The validity of our approach is tested on a conductance based neuronal model showing WLC and through the extraction of Lévy flights inducing fractals from recordings of rat hippocampus during open field foraging. Further insights are gained analyzing mice motor cortex neurons and non motor cell signals. The proposed mechanism provides a plausible explanation for the neuro-dynamical fundamentals of spatial searching patterns observed in animals (including humans) and illustrates an until now unknown way to encode information in neuronal temporal series.

The movement of many animals may follow Lévy patterns. The underlying generating neuronal dynamics of such a behavior is unknown. In this paper we show that a novel discovery of multifractality in winnerless competition (WLC) systems reveals a potential encoding mechanism that is translatable into two dimensional superdiffusive Lévy movements. The validity of our approach is tested on a conductance based neuronal model showing WLC and through the extraction of Lévy flights inducing fractals from recordings of rat hippocampus during open field foraging. Further insights are gained analyzing mice motor cortex neurons and non motor cell signals. The proposed mechanism provides a plausible explanation for the neuro-dynamical fundamentals of spatial searching patterns observed in animals (including humans) and illustrates an until now unknown way to encode information in neuronal temporal series.
In the context of animal movement the earliest reference to the superdiffusion of organisms is a study by Shlesinger and Klafter 1 about Lévy random walks (1986). Since this first theoretical prediction many studies have evidenced that the movement of several animals follow Lévy patterns. It is the case of movement patterns of wandering albatrosses 2 , jackals 3 , reindeers 4 , microzooplankton 5 , swarming bacteria 6 , spider monkeys 7,8 , root-feeding insects 9,10 , honey bees 11 , goats 12 , fresh water fishes 13 , snails 14 , fruit flies 15 , bony fish, sharks, sea turtles, penguins [16][17][18] , fallow deers 19 and humans 20 , just to name a few. According to optimal foraging theory 21,22 , evolution through natural selection has led over time to highly efficient searching strategies. Remarkably, it has been shown that for the case of revisitable targets Lévy search patterns optimizes the searching process 23 . Indeed, many works aim at identifying the ecological scenarios under which Lévy patterns emerge (or not) 18 . The study of search patterns in organisms is an active research area important not only in the field of movement ecology but also because of political, economic, environmental, and health-related reasons 24 . As a matter of fact, biological searching is a particular instance of a random search. Random searches can include quite different situations like searches performed by enzymes on DNA strands for specific sequences 25 . An authoritative account of this problem can be found in 24 .
A reduced number of works have advanced ideas regarding internal biological dynamics producing these search patterns. Examples include anomalous diffusion as a result of stochastic time delayed dynamics in the nervous system 26 , adaptive memory losses of previous behavior 27 or suppression of the scale-free power law behavior of fruit flies by blocking of synapses in the motor cortex 28 . In this context, one can ask if there are neuronal codification processes underlying the generation of animal search patterns. To the best of our knowledge this question is still unanswered, the underlying generating neuronal dynamics is currently unknown. In this work, we describe a plausible process for the neuro-dynamical origins of animal Lévy searches. The proposed mechanism is rooted in the so-called winnerless competition (WLC) dynamics, also known as cyclic dominance. We will show that WLC establishes a theoretical ground that facilitates information extraction from real neuronal dynamics that is translatable into Lévy search patterns.
Let us recall that WLC is grounded on heteroclinic cycles, i.e., on a collection of solution trajectories connecting equilibriums, periodic solutions, or chaotic invariant sets via saddle-sink connections. In WLC all participant agents alternate sequentially in time 29,30 in such a way that the system outcome can be considered a coding sequence 30 . WLC has been suggested as an archetypal dynamical principle useful to explain a diverse range of nervous responses [31][32][33] . WLC remarks the importance of itinerant and competitive dynamics in the brain as an information processing device. Experimental and theoretical work provides evidence that transient states can better represent information processing in the brain [34][35][36][37][38][39] . An interesting example of the application of the WLC principle is the case of the mollusc Clione limacina, where it is used to understand the competitive dynamics between statocyst receptor neurons and how these deliver control signals for swimming and hunting 40 .

Results
Detecting multifractality in WLC residence times. To set up the theoretical background on which this work is based let us consider the set of N coupled Lotka-Volterra (L-V) maps given by the equation where ρ ij represents the strength of the inhibitory action of i on j, > a 0 i represents the activity of each neuron, > r 0 is a control parameter and = , , …, i N 1 2 is the neuron index. Necessary and sufficient conditions for the existence of heteroclinic cycles supporting WLC for different values of N and r are known 41 . An illustrative temporal series showing WLC's characteristic sequential alternation of the activation between different neurons can be observed in Fig. 1a. To extract novel encoding possibilities from WLC let us proceed as follows: i) Given a particular neuron i and the time series (  , until an arbitrary visit index > n 0. An illustrative example of this series is depicted in Fig. 1b. ii) At this point, the differences between residence times of visits to the saddle, separated by a time interval , is calculated as   for arbitrary n. Note that these temporal series depend on the value of r. A particular result obtained applying scheme II is shown in Fig. 1c. The frequency content of the resulting oscillations is revealed by the power spectrum as is shown in Fig. 1d. This is a result obtained for a fixed value of the control parameter r. Therefore, we can calculate the power spectrum for all the values of the control parameter sustaining WLC behavior 41 , i.e., in the range < < r 1 3. The result of this calculation is displayed in Fig. 1e . It can be observed that for certain values of r there seems to be no frequency content, i.e., the differences between residence times are constant. It has been demonstrated 41 that residence times may show different regimens. In particular, they can display regimes of constant increments with the activation time (k), yielding constant differences between residence times. An additional situation producing constant differences corresponds to periodic residence times in which the constant equals to zero. Now, going back to Fig. 1e, it can be seen that Δ r-regions with nonzero frequency display bands containing Δ r-regions without frequency content. When observed in a smaller scale, a similar arrangement is found, evoking a certain ordering concealed on the ( , ) r f space (Fig. 1f). To analyze such a potential ordering we proceed by submitting the patterns on the ( , ) r f plane to a multifractal analysis. On the ( , )  Unveiling the neural code. We have been able to show hidden multifractality in WLC residence times. If transient sequential dynamics is subjacent to neuronal information processing 31 , it seems logical to us to enquire how this multifractality could encode neuronal information. We address this question in the following lines. Let us start at Fig. 2e, there we show the fractal dust, that we will denote by  , for a given > n 0. Such fractal dust has been extracted from the multifractal depicted at Fig. 2b, given an arbitrarily chosen fixed value of the control parameter, = .
1, carries useful functional information. With this idea in mind we can built a random walk on the plane such that Figure 2f shows the resulting random walk. It is characterized by a short scale spatial exploration alternated with jumps covering larger scales. These movement yields a Hurst exponent . > H 0 689 1 2 (Fig. 2g), i.e., it is a superdiffusive movement. We stress that not all sets ( ) r display the fractal structure required to yield a superdiffusive pattern. The reported pattern makes necessary tuning of the control parameter r.
As we said before, WLC dynamics is rooted in the existence of heteroclinic cycles. These trajectories can connect not only equilibrium points but also more complex solutions, as periodic or chaotic sets, producing temporal bursting behavior. We would like to know if the approach used above can also account for cases where bursting behavior is present. In such a case, instead of dealing with the difference between the residence times we would be dealing with the difference between the residence times in each one of the cycles composing a complex oscillation, i.e., the width of each constituent bursting spike. To better explain this case consider the temporal series obtained with Eq. (1) for the particular case of = .
r 2 962 (Fig. 3a). This signal shows evident bursting behavior. By establishing an arbitrary edge value (red line in Fig. 3a) we can determine the width of each constituent bursting spike and with this information we can proceed to calculate the difference between those widths (Fig. 3b) according to scheme II. After calculating the power spectrum for this temporal series (Fig. 3c) we extract a fractal dust, as shown in Fig. 3d. The fractal dimension of this dust is . d 0 873. Indeed, this fractal induces a bidimensional random walk (Fig. 3e) with a temporal evolution of its mean squared displacement characterized by a Hurst exponent . > H 0 858 1 2 and a distribution of steps characterized by an exponent µ − .
1 943 (Fig. 3g). This case exemplifies how a superdiffusive Lévy walk can be obtained from the underlying WLC multifractality. Note that the analysis of bursting spikes yields larger temporal series, therefore facilitating a proper estimation of the characteristic exponent. At this stage we have shown how a superdiffusive Lévy pattern can be obtained from WLC multifractality. However, our demonstration was not done using a realistic neuronal model. To overcome this limitation we consider a conductance based neuronal model also showing WLC dynamics 44 . Details about the model and its simulation can be found in the Supplementary Information. The conductance based model was analyzed with the same procedure as with the L-V map. Results for the scheme II are reported here. The temporal series obtained from the conductance model is shown in Fig. 4a and the series corresponding to the difference in residence times for visits separated a time interval = D 3, is depicted in Fig. 4b. The observed behavior is qualitatively similar to the one obtained from the differences in the case of the L-V map. The power spectrum for this signal is depicted in the Fig. 4c (calculated with 2 10 data points). We calculated power spectra for different values of the model parameter gaba and, following the same treatment as with the L-V model, the set of values ( , ) gaba f , satisfying ( ) = S f 0 gaba , are represented on the plane, as depicted in Fig. 4d. A multifractal analysis reveals multifractality for the set of points satisfying ( ) = S f 0 gaba (see Fig. 4e). In particular, for a fix value of gaba = 50 nsec, we can extract a fractal dust with dimension . d 0 612, as shown in Fig. 4f. This fractal allows for the construction of a random walk (Fig. 4g) whose mean squared displacement obeys a power law with a Hurst exponent , H 0 717. The distribution of steps does not show a promising shape. Even so, we fitted a power law to this collection of points obtaining an exponent µ − . 4 643. Clearly, this result is not conclusive of Lévy superdiffusion.
Decoding Lévy walks from experimental neuronal data. Notwithstanding these outcomes it is still unclear if real experimentally determined neural time series would share the same dynamics as the biomathematical models analyzed above. We address this concern analyzing temporal series from multichannel simultaneous recordings made from layer CA1 of the right dorsal hippocampus of three Long-Evans rats during open field tasks in which the animals chased randomly placed drops of water or pieces of Froot Loops while on a elevated square platform 45,46 . Here we show results from one of these recordings. The temporal series (see Fig. 5a) was processed following the scheme II. Figure 5 resumes our findings: it can be seen that the analyzed experimental data behaves as expected: The Δ -series displays complex oscillations (Fig. 5b) and exhibits a power spectrum (see Fig. 5c, it has been calculated with a temporal series of 2 17 points) whose zero power values allows us to extract a fractal set (Fig. 5d), with fractal dimension . d 0 707, that induces a Levy walk on the plane (Fig. 5e) (Fig. 5f) and a power law distribution of steps with an exponent µ − .
1 301 (Fig. 5g). Our approach seems to reveal an until now unknown coding mechanism based on a hidden multifractal. Indeed, it also points towards a relation between real neuronal dynamics and WLC, as evidenced by an excellent matching Scientific RepoRts | 5:18009 | DOI: 10.1038/srep18009 between theory and the outcomes obtained with the Long-Evans search related experimental temporal series. While these measures may not involve only motor related neurons we assume that the signal is carrying the relevant information related to the search process.
To extend our analysis we decided to further explore different available datasets. Consequently, we analyzed extracellular recordings from the anterior motor cortex neurons related to voluntary movement in mice 47 (see further details in the Supplementary Information). Results from this insight are summarized in Fig. 6. Again, a fractal set (this time with dimension . ) d 0 603 is obtained from the zero power values of the spectrum (calculated with a temporal series of 2 10 points), inducing a superdiffusive random walk with a Hurst exponent , H 0 702. When analyzing the steps' distribution it is found that the quality of the obtained profile is very poor. Even so, we attempt to fit a power law with exponent µ − .
2 708. Obviously, this outcome is far from conclusive of a Levy superdiffusion random walk. While this latter case is not fully satisfactory, it can be said that the present approach has been successful unveiling the fractal code hidden in the two motor neuron's series -i.e., both potentially related with searching -and obtaining (at least) superdiffusion. But, is this behavior exclusive of motor neurons or can it also be detected in non motor neurons -i.e., in neurons not directly related with searching -? To answer this question, we analyze available data from experiments with the grasshopper (Locusta migratoria) auditory receptor cell 48 (see Supplementary Information). Results for the grasshopper are shown in the Fig. 7 (the power spectrum was calculated with only 2 10 available points). The behavior of the mean squared displacement correspond to a superdiffusive regime characterized by a Hurst exponent , H 0 779 (Fig. 7f). However, the profile shown by the distribution of steps is insufficient to draw a solid conclusion. Again, we attempt to fit a power law, obtaining an exponent µ − .
2 498. This result reveals that the superdiffusive behavior is not exclusive to motor neurons recordings and that they can also be obtained from auditory receptor cells.

Discussion
It has been shown that WLC multifractality is an useful framework to generate neurologically based superdiffusive search patterns and to reproduce Lévy search patterns from a search related experimental temporal series (Long-Evans rats). Non conclusive results were obtained for Lévy search patterns in the cases of modeled and experimental neuronal time series not related with searching tasks. However, let's note that this result could be caused by the limited number of points available to generate a Lévy random pattern. The power spectrum's frequency precision is proportional to the number of data used in its calculation. In the cases of L-V and Long-Evans rats the temporal series were large enough to yield dense power spectra resulting in rich fractal dust. These sets were able to induce a random walk running on a diversity of steps lengths that were able to draw a well profiled power law. The detailed observation of the power spectra in Figs 4c,6c and 7c indicates that these power spectra do not show the frequency precision required to generate a fine grained fractal dust. Accordingly, one could expect to obtain Lévy search patterns also from larger temporal series in these situations. However, this can only be fully corroborated analyzing long temporal series. Thus, we have shown that Lévy search patterns can be obtained from non neurologically meaningful dynamics -the L-V equations-and real neuronal data related to a search process -the Long-Evans dataset-, while superdifussive random walks have been obtained in all the considered cases.
The conditions used in the Long-Evans (L-E) rat experiments do not seem to correspond to a search with revisitable targets as is required for an optimizing Lévy walk 23 . However, in this case once the rat drank the water the researchers dropped new drops randomly (György Buzsáki, private communication). While this fact does not exclude a particular spot to be filled again, chances are that this happened rarely. However, the animal has no information about this fact. Therefore, there is no reason to think that the animal is not revisiting small regions. It could happen if the water or food served in different instances fell randomly around the same area. Indeed, the rat would revisit such regions. Therefore we would not exclude revisiting dynamics in the L-E case. The biological mechanisms currently proposed to explain the generation of Lévy superdiffusion are not based on neuronal dynamics as the actual mechanism of WLC multifractality based coding does. The present approach may not show an evident relation with prioritization processes 49 or with bursts of rapidly occurring reorientations 50 ; two processes that have been related with the generation of superdiffusion 24 . In particular, in our case reorientation was dictated by an uniform distribution, i.e., unable to produce bursts of rapidly occurring reorientations. However, as it was already explained, in the current mechanism superdiffusion is generated by picking up randomly distances between the points of the fractal dust. This procedure may lead to a consecutive (bursting) extraction of large differences alternated with the consecutive (bursting) extraction of small ones. Therefore, one could conjecture the existence of some coincidences between the above mentioned ideas and the current mechanism.
In addition, a new approach to analyze coding in neurons has been presented. It may be possible to stimulate new experimental approaches to improve our insights into neuronal multifractal coding. The grasshopper results suggest that the underlying multifractality could be found at several processing levels, including receptor neurons. There is evidence that the metathoracic auditory system -to which auditory receptor neurons belong -could be a feedforward network 51 , however the presence of small but significant correlations between two receptor neurons does not allow to rule out the possibility of undetected weak synaptic contacts among receptors 51 , even though to present knowledge no synapses exist among them 52 . Our results favor the hypothesis that some sort of coupling between auditive receptors could be present. It seems important to recall pioneering results on the auditive receptor interaction of Locust that suggested electrotonic interaction among receptors in the region of the subreceptor plexus, given the absence of schwann sheaths and the penetration of receptor axons by collateral receptor axons 53 . Hence, it may be conjectured that auditive receptor neurons form a sensory network that support WLC multifractality. Notoriously, studies of sensory arrays mimicked by cellular automata excitable elements have shown collective nonlinear properties that improve input sensitivity and dynamic range 54 .
Instead of designing an experiment on the particular aspect of the generation of movement patterns perhaps it could be better to pay attention to the more general context of neuron coding. One may envisage that a possible experimental validation of this study could be achieved in a set up consisting of a neural system with a fractal dust precisely determined while the system controls a particular response. Subsequently, such a neural activity would be optogenetically turned off and the response triggered with a different simulated signal resembling the original one only in the fractal dust composition. An animal model for such a test could be the peripheral axons in the rat's sciatic nerve where optical inhibition of motor neurons and muscle activity in vivo 55,56 and in freely moving 57 non-transgenic animals has been demonstrated. A possible limitation with this idea is the requirement of real-time monitoring of the local field of a population of neurons in the targeted branch of the rat's sciatic nerve for its processing and fractal dust determination -this information would be used to produce the faked activation response after optical inhibition -. In this context it is still not clear how to achieve the process of faked signal injection or what would be the dimensionality of such a signal. No doubt the design of a proper experiment would need deeper insights to identify best fitted alternatives and/or possible simpler models. Summarizing, evidence of a WLC multifractality based coding mechanism was presented. This phenomenon unveils a parameter-tuned fractal set that induces a 2D superdiffusive Lévy random walk for some of the cases analyzed. The current insight provides a theoretical ground to analyze modeled and real experimental neuronal time series. This work shows evidence about the possible neuronal basis of observed Lévy search patterns in animals. It should be remarked that the observed patterns are a result of tuning a control parameter. There are control parameter values for which a fractal dust able to induce superdiffusion can't be obtained. We call attention to this point as it may be related to the fact that certain animals may not be following superdiffusive Lévy search patterns 18,24 or displaying a multi-scale pattern 58 .
One could conjeture that if Lévy search patterns are an optimal searching solution, it could also be present in the computational exploration of a space of solutions associated to complex nervous system tasks. Obviously, a test of this conjeture requires further and deeper insights. Finally, it is well known that each spike carries information 59 . Remarkably, our approach is able to extract the information contained in the width of each spike -the duration of each bursting spike-to translate it into a detectable code. In this context, it must be noted that detectability is related to the determination of the constituent frequencies. However, from the point of view of nervous information processing, operating conditions based on long temporal series may be impractical. This suggest that the encoding mechanisms presented here may be well suited for better performance in a network context were each neuron generates its own (low-intermediate precision) fractal dust based on short temporal series, while the full network could induce well profiled Lévy walks superimposing all the generated single fractal dust. Our understanding of the discovered mechanism is still in an initial stage implying that it may contain a certain degree of speculation. Maybe the observed multifractality is also present in different nonlinear dynamical scenarios with the same coding potential.

Methods
Simulation of the L-V map was done using the connectivity matrix specified in the Supplementary Information. The conductance model was numerically integrated with a backward differential formula (more details in the Supplementary Information). Peaks widths for the models and the experimental time series were determined fixing the threshold maximizing the number of widths. The Hurst exponents were determined using the Generalized Hurst exponent approach 60 and the mean squared displacement was plotted based on the displayed random walk. Figure 6. Superdiffusive random walks decoded from mice motor neurons. We extend our analysis to include the (a) action potential from adult mice (used threshold in red) whose (b) ∆ ( = ) D 2 i signal also behaves in a complex manner displaying a (c) power spectrum calculated with 2 10 points with a diverse frequency content. From values of frequency satisfying ( ) = S f 0 we can determine (d) a fractal dust whose gaps induces a (e) 2D random walk with (f) a superdiffusive mean square displacement calculated as in Figure 2(g). (g) The number of available data points used in the power spectra calculation is not enough to detect the possible Lévy character of this superdiffusion (best fit in red).