Statistics of Decision Making in the Leech

Animals and preparations

Adult leeches (Hirudo medicinalis) were obtained from Ricarimpex (Eysines, France) and kept at 5°C in tap water dechlorinated by aeration for 24 h. Two different kinds of preparations were used. The first preparation consisted of an intact leech free to move on a dish. The second preparation consisted of an intact leech restrained by fixing the rear sucker of the animal to the center of the dish with one pin. Pins were carefully inserted through the animal's body, avoiding the penetration of the connective fibers running along the animal. When one pin was used, the leech was able to fully rotate around the pin and could explore the surrounding space.

All preparations were kept either in a round dish (15 cm diameter) or a square tank (15 × 15 cm) filled with dechlorinated water at constant temperature (20-24°C) and illuminated by dim, diffused light, with no abrupt spatial and/or temporal gradients of sensory stimuli, such as chemical or olfactory anisotropies, isolated visual cues, or abrupt illumination changes. The dechlorinated water in the dish was changed every 2 h to avoid the effects of serotonin and ions released by the leech (Schnizler and Clauss, 1998; Schnizler et al., 2002). In every preparation, three colored plastic beads (red, green, and blue) of a diameter of ∼2 mm were glued on the back of the leech with Superglue (World Precision Instruments, Berlin, Germany). The red, green, and blue beads were attached near the head, the middle, and the tail of the leech.

Imaging

A color CCD camera (640 × 480 pixels of image size; model 231S; Watec, Tsuruoka, Japan) was used to view the leech from above. The camera was connected via the S-video-output to a frame grabber (PCI-1411; National Instruments, Milan, Italy) installed on a personal computer, which processed the images in real time. The colored beads placed on the back of the leech were tracked at 10 Hz by using a software program developed by Sergio Graziosi in our laboratory (Mazzoni et al., 2005) with Lab-VIEW 6.1 (National Instruments). Briefly, the software acquires images directly in the hue/saturation/lightness color space, and the tracking algorithm works in real time with no frequency constraints. A user-customized color matching is performed to convert each acquired frame into three binary images, one for each bead, in which the ON bit represents a pixel that fit in the expected color subspace. This software did not require saving all of the processed images but only the coordinates of the selected beads. In this way, it was possible to monitor the leech behavior for several hours without filling hard disk space. The positions of the beads on the image plane were acquired as Cartesian coordinates, and six time series were obtained as follows:

where every step corresponds to 100 ms.

Some analysis required the motion to be studied in a polar system centered on the middle of the dish. In this way, Cartesian coordinates (x,y)_color were transformed into a pair of planar polar coordinates (ρ,θ)_color. When the leech was restrained from moving around by a pin inserted in its tail, the polar system had a clear geometrical and physiological meaning: ρ_color is the elongation of the leech body at the position marked by the bead (color), and θ_color is the angle between this position and a reference axis.

Sensory stimulation

In one series of experiments, a small pulsating light-emitting diode (LED) (at 0.5 Hz) was put above the water but near the edge of the dish. The behavior of the leech was first monitored for 3 h in control conditions of low and diffused light (without light stimulation). The animal was allowed to rest for 30 min, and then the light stimulation (LED) was initiated. In some experiments, two LEDs were placed diametrically opposite each other and flashed alternately at 1 s intervals. These experiments were done with restrained leeches, and the residence time distributions, in both light and control (no light) preparations, were analyzed.

Behavior classification

Often the leech moved in the dish for >12 h. The marker coordinate time series thus often comprised >100,000 samples, which had to be analyzed consistently. We therefore developed an automatic classification method for identifying the leech behavior, which is based on analyzing the kinematics and spectral content of the recorded time series and is described in detail in a previous work (Mazzoni et al., 2005).

Briefly, leech movements were classified into six different categories (Table 1): stationary (including still, peristaltic-like, and head-attached states), swimming, pseudo-swimming (or ventilation states), crawling, and exploratory and abrupt movements. At each time step n of the recorded time series, different criteria were used to discriminate one behavioral category from the others. For example, given the set of coordinates describing the leech dynamics at time step n, to test whether the current behavioral state was a stationary state, the speed of the head [V_head(n)] and tail [V_tail(n)] beads at time n was calculated. If the speed was less than 1 pixel per second (pixel/s), then the state was considered stationary. Additionally, to distinguish between the three different stationary states, the speed of the midbody [V_midbody(n)] bead and its frequency were computed. If V_midbody(n) was <1 pixel/s, then the state was classified as still; however, if V_midbody(n) was >1 pixel/s, then the frequency of oscillation was computed, and, depending on this result, the stationary state was classified as peristaltic like (with a frequency of ∼0.03 Hz) or head attached, during which the midbody oscillated in an irregular way.

The other behavioral categories were identified by computing the following parameters: (1) elongation e(n), that is, the distance from head to tail beads; (2) maximal speed of the head over a time window of width T centered in n (with T ranging from 20 to 50 s); (3) dominant frequencies of the power spectrum of e(n), computed over the same time window of width T. More specifically, oscillatory behaviors (such as swimming, pseudo-swimming, and crawling) could be reliably identified by the pairs of values [, ]. For instance, swimming episodes were consistently characterized by and , whereas crawling was characterized by and . The spectral analysis of e(n) also allowed identification of the exploratory states, because events with at the lowest resolved frequencies (∼0.1 Hz) matched correctly with the exploratory movements (i.e., slow and irregular oscillations of the head and the anterior part of the body). An exhaustive inventory of all identified leech behaviors and the criteria used to classify them is given in Table 1. Note that all definitions of behavior are based on units of pixels, in which 1 pixel = 0.4 mm.

Statistical analysis of stationary, swimming, and exploratory states Space distribution. Once a specific behavior was identified across the whole time series, the mean values of the polar coordinates for the head [ρ_behav(i), θ_behav(i)] in each time interval i = 1,..., N_behav in which the leech was performing that behavior were computed (N_behav is the total number of episodes found for that behavior). The distributions of angles θ_behav(i) were plotted as polar histograms (see Fig. 2 B, D), with an angular bin size of 18°. If the analyzed behavior is isotropic, then the number of recorded states with θ_behav(i) belonging to a given angular sector must be distributed according to a Poisson distribution with mean μ = (number of occurrences N_behav)/(number of sectors S). A χ² test was used to check whether the obtained distributions were Poisson distributions (p > 0.05) and therefore whether the behavior was isotropic. If the behavior was found to be significantly anisotropic, then a post hoc test was used to specifically test whether such anisotropy was caused by one or more sectors that were overpopulated and to find out the orientation of such sectors (i.e., the preferred direction for that behavior). The post hoc test was based on the null hypothesis that the distribution of states per sector is not significantly different from a Poisson distribution and that therefore the probability of having Z states or less in a sector is the following:

Then, the probability Q that n exceeds Z, at least in one sector, is the following:

In Figure 2 B-D, the magenta circle indicates the value of Z associated with Q = 0.05. Note that this value depends only on the number of states and the number of sectors.

Time series correlation. To study the linear properties of the sequence of identified states, series of average angles θ_behav(i) and lengths ρ_behav(i), with i = 1,..., N_behav, were considered as the outcome of two stochastic processes Θ_behav(i) and P_behav(i). First, we assessed whether these stochastic processes were stationary, testing the stability of the mean and variance on the first and second half of the series (ANOVA test).

The correlation between consecutive values θ_behav(i) [and consecutive values ρ_behav(i)] was then evaluated by computing the autocovariance of the process Θ_behav(i) [and of the process P_behav(i)]. Similarly, the correlation between consecutive pairs of values θ_behav(i) and ρ_behav(i) was evaluated by computing the covariance between the processes Θ_behav(i) and P_behav(i). Statistical significances of autocorrelation and cross-correlation values were evaluated by a two-tailed normal test under the assumption that correlation values for uncorrelated time series are normally distributed with mean zero and variance 1/(N - l - 3), where N is the number of samples in the two series, and l is the lag number. Correlation values were considered significant at p < 0.05. The number of consecutive significantly correlated states around n = 0 gives the size of the central peak of the autocorrelation. The width of the central peak of the autocovariance function gives a measure of the “memory” of the stochastic process and indicates on how many previous states the identity of the state depends (Priestley, 1981). To compare the processes Θ_behav(i) and P_behav(i) with random processes of the same statistical properties but without any correlation (“no memory”), the original series ρ_behav(i) and θ_behav(i) were scrambled before the autocovariance was computed.

Time structure of the sequence of behaviors. The time structure of the sequence of behavioral states was studied by computing, for any identified behavior, the distribution of the duration Δt of the time intervals (residence times) during which that behavior was performed. In every experiment and for every behavior, this distribution had a high and narrow peak centered on a small Δt (usually <10 s), which decayed to 0 with increasing Δt. This suggested that the sequence of behaviors could be modeled as a Markov process in which the residence time for each behavior is distributed as an exponential or a mixture of exponentials (Haccou and Meelis, 1992). Therefore, for each behavior, the observed residence time distribution was fit by an exponential function, f(Δt) = Aexp(-Δt/τ), or by the sum of two exponential functions. To determine whether two exponentials were required, goodness of fit was tested with a χ² test, and a one-exponential fit was accepted only if there was no significant difference (p > 0.05) between the data and the fit. When one exponential was insufficient, two exponentials were always found to be sufficient by the χ² criterion.

Power spectral analysis. An informative parameter about the motion of the leech is the elongation e(n), i.e., the distance between the head and the tail, which varies during every oscillatory behavior. As mentioned previously, during the classification procedure (Mazzoni et al., 2005), for every nonstationary point in the time series, a power spectrum PS(f)of the elongation over the portion of the time series centered on that point and having a half-size of 10 s was computed to accurately detect high-frequency oscillations. A second power spectrum with a half-size of 25 s was then computed to discriminate low-frequency oscillations. Sometimes, power spectra were computed even for much larger windows to analyze global properties of the behavior (see Fig. 4 A, B). Welch's averaged modified periodogram method was used (pwelch function in Mat-lab; MathWorks, Natick, MA), with 50% overlap between successive time windows, each one being one-fourth the size of the analyzed interval. Points of the resulting power spectra were averaged within logarithmically spaced frequency bins to smooth the power spectrum at high frequencies. Power spectra were further smoothed by averaging every spectral point PS(f) over a semi-octave centered on each frequency f to enhance detection of any peaks (Mazzoni et al., 2005).

Recurrence analysis of exploratory states. Complex exploratory motion was investigated using the recurrence plot technique (Eckmann et al., 1987). The advantage of this approach is its capability of detecting deterministic or repeated component patterns even in highly nonlinear and irregular time series, which would be lost through averaging in linear measures such as the power spectrum (or equivalently, the autocorrelation function). First, an m-dimensional delay representation (Kantz and Schreiber, 1997) is constructed for a 1000 s segment of motion time series. The delay coordinate vectors d_n are given by the following:

These represent successive short pieces of the waveform history preceding each point y_n in the time series. We used m = 20 and lag v = 2. A 4000 × 4000 diagonally symmetric matrix c is then plotted in which the point at row i, column j is color coded (in gray scale) according to the Euclidean distance between d_i and d_j as follows:

Within this plot (see Fig. 5), dark colors represent similarity, or low distance, between delay vectors. A sequence of intensely dark points along a line of slope -1 represents a flow in which history repeats itself. We visually scanned the recurrence matrix for diagonal lines satisfying a threshold similarity condition of length at least twice the duration of the delay vector (>19 s). Having identified one diagonal, all similar repeats within the time series are easily extracted as those diagonal lines that intersect a horizontal (or vertical) line drawn through the center of the first line. The repetition frequency of the pattern is quantified by the number of such sister diagonals (for a diagram of this process, see Fig. 5). Repeats that follow the same trajectory but with a different timescale appear as diagonal lines with different slopes.

Estimation of transition probabilities and test of sequential dependences. Once every point of the time series has been classified, the behavior of the leech can be seen as a sequence of identified behavioral states. To estimate the transition probabilities between pairs and triplets of behaviors, we computed, for each behavior A, the following quantities: the number of occurrences N_A of behavior A; the number of times N_AB that behavior A was followed by behavior B (i.e., the number of transitions A → B), for all possible following behaviors B; the number of times N_XA that behavior A was preceded by behavior X (i.e., the number of transitions X → A), for all possible preceding behaviors X; the number N_XAB of transitions X → A → B, for all possible preceding behaviors X and following behaviors B. Then, the probabilities of first-order transitions A → B were estimated as P(B/A) = N_AB/N_A. Similarly, the probabilities of second-order transitions X → A → B were estimated as P(B/A, X) = N_XAB/N_XA. Values of these transition probabilities during the first and second half of the experiments were compared with a χ² test to check whether the process was stationary. Differences between the two sets were considered significant at p < 0.05.

To test whether the sequence of behavioral states could be modeled as a first-order Markov chain, we used the χ² test for first-against second-order dependences in the sequences of identified states (Haccou and Meelis, 1992). Given the null hypothesis that P(B/A) = P(B/A, X), for all triplets X → A → B, we first computed the χ² test statistics for each behavior A: , where the summation is taken over all behaviors X and B, preceding and following A. Then, a global χ² test statistics was computed as χ² = ∑_A χ²_A, i.e., taking into account all possible intermediate behaviors A. The degrees of freedom of the global χ² test statistics are k = b(b - 2)², where b is the number of different behaviors identified in the analyzed experiments; however, if some of the transitions toward or from A occur rarely or not at all (i.e., with almost zero probability), this number is overestimated. In this case, the correct number of degrees of freedom becomes k =∑_A(b - v_A - 1) (b - w_A - 1), where v_A and w_A are the number of transitions toward A and from A, respectively, that cannot occur (Haccou and Meelis, 1992). Deviations from the null hypothesis of a first-order-dependent Markov chain were considered significant at p < 0.05.

Behavior inventory and distribution of the durations of each state

Figure 1A illustrates the time evolution of the elongation e(n), i.e., distance between the head and tail bead, while the leech moved freely in a covered dish. The analysis of e(n) (Mazzoni et al., 2005) allowed us to identify well known leech behaviors such as crawling and swimming, which are shown in Figure 1A (example of swimming, dark gray, and three superimposed snapshots, right).

By analyzing the motion of the head, midbody, and tail, some other behavioral categories were identified (Mazzoni et al., 2005) and are summarized in Materials and Methods and Table 1. Additionally, there were some unclassified states, characterized by slow drifts or by oscillations superimposed on an apparently exploratory event. One example of this is shown between two swimming states in Figure 1A (white). Unclassified states comprised, on average, ∼10% of the observation time.

Leeches not only show stationary and periodic states, but they also exhibit a more complex behavior, referred to here as exploratory (Fig. 1A, light gray). This behavioral category is characterized primarily by the absence of any evident regularity (Table 1). Three superimposed snapshots of exploratory behavior are shown in the left inset of Figure 1A. Figure 1B shows the relative occurrence of different behavioral categories observed over 6 h of continuous recording of leech motion (three restrained leeches). For each leech, statistics were obtained from 0 to 3 h (left bar) and from 3 to 6 h (right bar), respectively. Figure 1C shows the same statistics but for three freely moving leeches. These data indicate a degree of variability among individuals, but the repertoire of behaviors identified for a given animal as well as the percentage of time spent in each behavior were well preserved during the first and second 3 h of observation.

The distribution of the residence times of exploratory, swimming, stationary, and crawling states is shown in Figure 1D-G. The distributions for exploratory and swimming states could be reliably fitted (χ² test; p > 0.05) by a single-exponential function with time constants τ = 26.8 s and τ = 42.5 s, respectively (Fig. 1D,E). On the other hand, the distributions of the residence time for stationary and crawling states were significantly different (χ² test; p < 0.05) from the best fit achieved with a single-exponential function; however, in both cases, the sum of two exponential functions gave a best fit that did not significantly differ from the observed distribution (χ² test; p > 0.05). The time constants found were τ₁ = 4.5 s and τ₂ = 40.5 s for stationary states and τ₁ = 5 s and τ₂ = 58 s for crawling states (Fig. 1F,G). All of these fits have been calculated on distributions of pooled states from many leeches (7 for exploratory, 4 for swimming, 11 for stationary, and 4 for crawling) because single experiments did not have enough data to provide a reliable evaluation of the time constants; however, after the best fit for each pooled distribution was found, we tested whether individual distributions were consistent with that fit. The results of this second test are as follows: six of seven residence time distributions of exploratory states, four of four distributions of swimming states, 11 of 11 of stationary states, and four of four distributions of crawling states were not significantly different (χ² test; p > 0.05) from the corresponding fitting functions obtained with the pooled data.

Spatial distribution

Leeches appear to wander around uniformly. Indeed, when the tail is pinned to the center of a round or square dish, the leech explores uniformly in every direction (Fig. 2A), and the orientation of its body during stationary states is evenly distributed over all θ coordinates. To test the isotropy of leech orientation during stationary states, the area of the dishes was divided into 20 sectors of 18° each (Fig. 2B). To have expected frequencies of at least five states per sector, seven experiments were done with pinned leeches, with at least 100 stationary states each, were analyzed. In seven of seven leeches, the distribution of the number of states per sector was not significantly different from a Poisson distribution (χ² test; p > 0.05), as expected for an isotropic behavior. As an example, four superimposed representative polar plot distributions (a different color for each leech) are shown in Figure 2B. These distributions show the leech orientation for the first 100 stationary states. The magenta circle indicates the confidence threshold for each sector, with p < 0.05 of being exceeded (see Materials and Methods), given the null hypothesis of isotropic behavior. The radius of this circle depends only on the number of states analyzed and on the number of sectors and therefore is the same for every data set analyzed in the figure. If one sector exceeds this confidence threshold, then the distribution of states is anisotropic, because the animal visited that sector a significantly higher number of times than the others. For the experiments shown in Figure 2, because no sector exceeded this threshold, there was no significant deviation from isotropy. Similar distributions were obtained for all experiments (n = 7).

Free leeches spent most of the time near the dish edges (Fig. 2C), but again the θ distribution of stationary states did not have any statistically privileged direction. For this analysis, five experiments with free leeches with at least 80 stationary states were analyzed (just three of them having >100 states). The number of occurrences per sector followed a Poisson distribution in five of five leeches (p > 0.05). As an example, the superimposed polar distributions for the first 50 stationary states observed in four different leeches (to simplify the plot) are shown in Figure 2D with their confidence threshold.

Temporal distribution

Elongations ρ and orientations θ during stationary states were ordered according to their temporal occurrence, and their temporal statistics were evaluated. Autocorrelation analysis for both ρ and θ (Fig. 2E,F, black line) exhibited a small but significant bump (normal test; p < 0.01) flanking the central peak, with a width corresponding to two to eight successive states that were significantly correlated. This bump was not observed when the original time series was scrambled (Fig. 2E,F, gray line). This phenomenon was observed and was significant in seven of seven free leeches and in seven of seven pinned leeches, with an average central peak width for ρ series (Fig. 2G, black bars) of 5 ± 2 states for pinned leeches and 4 ± 2 states for free leeches. The same results were obtained for θ series (Fig. 2G, white bars), with a central width of 6 ± 2 and 3 ± 1 states, respectively, indicating the existence of a correlation between positions of successive stationary states. Values of ρ and θ were not significantly cross-correlated.

Effects of sensory stimulation

We further tested how the spatial distribution of leech locations changed as a consequence of localized sensory stimulation. For this purpose, we used seven pinned leeches (different leeches from those in Fig. 2). Five of them were monitored when one light was flickering, and the other two were monitored when two lights were flickering. Figure 3 shows some representative polar plots in which the angular distribution of the position of the leech during stationary states is reported. Without sensory stimulation, the distribution of stationary states for all tested leeches was isotropic (the number of occurrences per sector followed a Poisson distribution; p > 0.05). As an example, the distribution of stationary states for three leeches is shown in Figure 3, A, C, and E, in which no sector exceeds the significance threshold (dark gray circle) for anisotropy; however, in four of five cases with one light and in two of two cases with two lights, the presence of the light broke this isotropy. Figure 3, B, D, and F, shows the angular distributions during stationary states for these same three leeches when the lights were switched on (the gray triangles indicates the orientation of the lights). In all three cases, some sectors were more populated than the statistical threshold for acceptance of isotropy (p < 0.05) (dark gray circle). Although a systematic relationship between orientation of light and preferred resting orientation of the leeches could not be found, in general the leeches avoided spending time under the light sources (Fig. 3D). Their preferred direction was often close to but not coincident with the light source, as if leeches moving toward the light source were constrained to proceed along a nearby direction (Fig. 3B).

Analysis of exploratory states

Well resolved peaks characterized the smoothed power spectrum of swimming and crawling episodes, at 1.5 and 0.2 Hz, respectively (Mazzoni et al., 2005). On the other hand, exploratory states were characterized by a 1/f^S distributed power spectrum (where S is the slope of the power spectrum when plotted on a logarithmic scale), with no resolvable peaks. Furthermore, the shape of the smoothed power spectrum of e(n) during exploratory episodes was very similar across different free and pinned leeches. Examples showing a superposition of four power spectra corresponding to different exploratory episodes observed in four pinned (reddish curves) and four free (bluish curves) leeches are reported in Figure 4, A and B, respectively. Similarity among exploratory power spectra was better quantified by computing their 1/f^S fit (Fig. 4A,B, black lines). The slope S of the power spectrum was very similar across different pinned (red bars) and free (blue bars) leeches, as shown by the slope distributions reported in Figure 4C. The value of the slope S is highly indicative of the kind of process analyzed and can be considered, in this case, as a measure of the correlation of the different movements composing the behavior (Berg, 1993; Szendro et al., 2001). The observed similarity suggests that, although very different and very irregular, all exploratory episodes belong to the same dynamics. For instance, the average slope of the power spectrum was 1.8 ± 0.2 for free (blue bars) and 2.0 ± 0.2 for pinned (red bars) leeches. Although small, the difference between average slopes in the two sets was significant (ANOVA test; p < 0.05). The slope of the power spectrum is close to 1/f², indicating that exploratory movements approximately follow random walk dynamics (Berg, 1993).

Despite the quasi-periodic appearance of exploratory motion (Fig. 1A, first light gray patch), its power spectrum was not associated with peaks at any frequency. Such a broadband spectrum is indicative of noise or nonlinear deterministic chaos. As a test for low-dimensional nonlinear deterministic chaos, which could be caused by, for example, the interactions of a few underlying motor oscillators, we examined whether exploratory time series could be predicted better by linear prediction, by nonlinear prediction, or by nonlinear prediction of the phase-randomized time series (Kantz and Schreiber, 1997); however, the best result was obtained using linear prediction (data not shown). Thus, we were not able to find evidence for simple low-dimensional chaos in exploratory time series.

We also verified that, during exploratory events, the leech visited the entire surrounding environment (Fig. 4D), with no preferential direction or location, and that the distribution of the head speed V_head was Gaussian (Fig. 4E), as it was for the other behaviors (Mazzoni et al., 2005). Even if there was a characteristic speed for this behavior (μ = 35 pixels per second) that was similar for pinned and free leeches, the distribution was broad (σ = 13 pixels per second) compared with other behaviors (Mazzoni et al., 2005), as shown by the statistics for seven pinned and seven free leeches presented in Figure 4F. This broad distribution suggested heterogeneity of activated motor patterns, as expected for the highly irregular motion characterizing the exploratory behavior. Indeed, using the recurrence plot technique, it was possible to identify qualitatively various repeated, effectively deterministic motifs within the exploratory motion. Figure 5A shows an example of a recurrence plot of a 1000 s section of exploratory behavior. A pattern of similar behavior repeated at two different points in the time series was observed as a diagonal line segment of slope -1. All other similar patterns within the time series could be picked out as those diagonal line segments, which intersect the same horizontal as the first identified diagonal. An example of this is shown in the inset; this pattern was found several times along the recurrent plot, and some of these are shown superimposed in Figure 5C as pattern 4. This pattern and some other recurrent patterns (Fig. 5C, patterns 1, 11) are superimposed in black on a piece of the original time series in Figure 5B. These patterns show that the evolution of the exploratory behavior is similar at different times over periods lasting tens of seconds.

We also performed the same analysis after Fourier transforming, randomizing the phase information, and taking the inverse transform. This phase-randomized surrogate series (Kantz and Schreiber, 1997) has the same power spectrum and autocorrelation as the original time series and therefore should give the same results if a linear stochastic process produced the data; however, there were in the original time series a number of recurrent patterns that were not observed in the phase-randomized surrogate. The 12 most common of these are shown in Figure 5C. They could be separated into several groups: (1) abrupt movements, as patterns 1 and 4, indicated in the original trace shown in Figure 5B, and (2) quasi-periodic fragments and oscillating extensions sandwiched by smooth oscillations, such as patterns 7-12.

Transitions between behaviors

A process in which the probability of transition between different states is determined only by the identity of the most recent state is said to be markovian. The exponential distribution of duration of exploratory, swimming, crawling, and stationary states (Fig. 1D,G) shows that every behavior or subcategory of behavior had a probability of ending per unit of time equal to 1/τ (where τ is its characteristic time constant). This is highly suggestive of a Markov process underlying the decision mechanisms that control transitions among different behaviors. To verify this possibility, transitions between different behavioral categories were analyzed. The probability P(B/A) of the first-order (or one step) transition A → B, i.e., the probability of state k + 1 belonging to behavior B if state k belongs to A, was estimated for each pair of behaviors (A, B), as described in Materials and Methods. The corresponding transition rate, i.e., the probability per unit of time for the leech of switching to behavior B when it is in a state belonging to the behavior A, is then given by P(B/A)/τ_A.

Figure 6A shows, for a pinned leech, one-step transition probabilities between five different behaviors, for the first (0-3 h) and second (3-6 h) half of the experiment (black and white bars, respectively). In each panel, the transition probabilities in the first half of the recording session were not significantly different from those of the second half (χ² test for homogeneity of distributions; p > 0.05), and the same result was obtained for four of five leeches examined. This indicates that the decision process was stationary in time.

Transition probabilities varied from leech to leech, but some properties of the Markov process were present in all leeches. For instance, peristalsis was always followed by stationary states (Fig. 6A, fourth panel), and transitions usually occurred among stationary, exploratory, and stationary states (Fig. 6A, first three panels). Some transitions, like swimming to crawling and crawling to swimming, were never observed in free leeches (six experiments); however, when the water level was not sufficient to fully cover the leech body, hybrid behaviors such as swimming-like oscillations superimposed on the elongation phase of crawling (Fig. 7) (Esch et al., 2002) and elongation followed by a swimming behavior were sometimes observed (data not shown). Figure 7A shows the time displacement of the x and y coordinates of the head (red), midbody (green), and tail (blue). Four crawling steps are shown during which the body was oscillating. In Figure 7B, the space displacement of the same hybrid behavior is shown. As can be seen, the leech was crawling from the middle to the bottom of the plot (y direction) with oscillatory movements (x direction) superimposed on the elongation step. Note the difference between the x and y scales.

To verify that the sequence of behavioral states could be modeled as a first-order Markov process, we tested whether transitions between different states satisfied the Markov first-order dependence conditions, i.e., P(B/A,X) = P(B/A), for every triplet X → A → B, where P(B/A,X) is the probability of state k + 1 belonging to behavior B if the current state k belongs to A and the previous state k - 1 belonged to X.

A χ² test (see Materials and Methods) was used to test whether the set of two-step transition probabilities P(B/A,X) was significantly different from the set of one-step transition probabilities P(B/A). The χ² statistics were computed taking into account all possible transitions between the three most frequent behaviors (stationary states, exploratory, and swimming) (Fig. 6A) in the five experiments that had an appreciable number of transitions (n > 5) expected for every triplet of states (for instance, Sw → Exp → Stat or Exp → Sw → Exp, in which Sw is swimming, Exp is exploratory, and Stat is stationary). The result of the test was that, for five of five experiments, the two-step probabilities were significantly different (p < 0.05) from the one-step probabilities, and therefore the sequences of identified behaviors could not be considered as first-order Markov processes; however, a comparison among the different terms in the χ² statistics revealed that deviations from first-order dependence were attributable mainly to the transitions Exp → Stat → Exp (Fig. 6B, first panel). This finding, together with the previous result that the distribution of the residence time for stationary states could be fit by a sum of two exponentials (Fig. 1F), strongly suggested that states identified as stationary were actually divided into two distinct subcategories. Therefore, based on the residence time distribution, we divided stationary states into two subcategories: short-lived stationary (SLS) states (duration <10 s) (Fig. 1F, first data point) and LLS states (duration >10 s).

The χ² test was then repeated to include those transitions involving only SLS states or those involving LLS states. When only short stationary states were included, the two-step probabilities were still significantly different (χ² test; p < 0.05; five of five experiments) from one-step probabilities. This again was because P(Exp/SLS, Exp) was always higher than P(Exp/SLS) (Fig. 6B, second panel). On the other hand, when only long stationary states were included in the statistics, conditioned one-step and two-step probabilities were not significantly different in four of five experiments (χ² test; p > 0.05). A comparison between P(Exp/LLS, Exp) and P(Exp/LLS) is shown in Figure 6B (third panel).

Overall, this analysis suggested that LLS and SLS states had very different properties. Long stationary states behaved like exploratory and swimming behaviors, and the transition probabilities were not affected by the previous behavior. On the contrary, during a short stationary state, the leech remembered the previous behavior and had the tendency to repeat it. This happened mainly when the leech was exploring. In fact, many of the short pauses were between two exploratory states, and most of them adjoin at least one exploratory state (32 and 74%, respectively, of SLS states for nine leeches in which both SLS and exploratory states had an appreciable number of occurrences). This suggests that SLS states are pauses that do not have to be considered as independent from the exploratory states. Therefore, we repeated the statistical test for first-order Markov dependence by absorbing each SLS state that bordered an exploratory state into the exploratory state itself. Such paired states contributed, on average, 15% of the total number of stationary states and 24% of the number of exploratory states (nine leeches). Figure 6C shows some of the resulting two-step and one-step transition probabilities after this adjustment in the classification of behaviors. Here, transitions from exploratory (Fig. 6C, top) and from swimming (bottom) to LLS states are compared, when the state preceding the swimming or exploratory state is either arbitrary (black) or LLS (gray). For four of five experiments, two-step probabilities were not significantly different (p > 0.05) from the one-step probability set, indicating that the underlying transition dynamics behaves, in essence, as a first-order Markov process.

Flow of leech behavioral decisions as a first-order Markov process

The results shown in this study (Figs. 1, 6) suggest that the time structure of leech behavior is well described as a first-order Markov process with transitions among various different behavioral categories, as illustrated in Figure 8A. Transition probabilities among different behavioral categories are likely to depend on the spontaneous neural activity of the animal brain as well on proprioceptive and sensory information. The discovery of many command interneurons in the subesophageal ganglion of the leech, the locations of which are indicated in the diagram in Figure 8B, offers a neurobiological basis for the scheme of Figure 8A.

Neurons Tr1, SE1, and R3b1 initiate swimming (Brodfuehrer and Friesen, 1986a,b,c; Brodfuehrer et al., 1995; Thorogood and Brodfuehrer 1995; Shaw and Kristan, 1997; Esch et al., 2002), whereas neurons SIN1 and Tr2 terminate swimming (Brodfuehrer and Burns, 1995; O'Gara and Friesen, 1995). The firing of neuron R3b1 (Esch et al., 2002) elicits both swimming and crawling or a hybrid behavior (similar to that shown in Fig. 7), depending on the sensory environment (the depth of the fluid surrounding the leech); however, leech behavioral decisions are likely to depend on a much more distributed neural dynamics than just the activity of this set of neurons. In fact, behavioral choice experiments have shown that most of these neurons are multifunctional, i.e., they are activated during different behaviors such as swimming and whole-body shortening (Tr1 and SE1) (Shaw and Kristan, 1997) or swimming and crawling (R3b1) (Esch et al., 2002). Moreover, recordings of the neural activity descending from the head ganglion in the connective have shown a massive coactivation of many unidentified interneurons during swimming episodes (Brodfuehrer and Burns, 1995). Finally, direct stimulation of identified command neurons produces highly variable behavioral responses and often does not elicit any kind of behavior (Brodfuehrer and Friesen, 1986b). The same variability is encountered when strong mechanosensory stimuli are used to elicit swimming (Grobstein, 1994; Brodfuehrer and Burns, 1995). As a consequence, the resulting stimulus-response curves do not have a clear stimulus level threshold for reproducible initiation of swimming. Overall, these findings suggest that leech behavioral decisions are controlled by the dynamics of a distributed neural network that defines the internal state of the leech CNS. This dynamics, because of its distributed nature and the intrinsic variability of leech interneuron and motoneuron firing (Zoccolan et al., 2002), can be regarded as a temporally structured probabilistic process. Therefore, we believe that the approach presented here (i.e., modeling leech behavioral data as a first-order Markov process) is at present the most effective for quantifying and characterizing the complexity of decision making in the leech CNS.

A future challenge in the study of leech behavioral decisions will be to compare transition probabilities in leeches moving in an unstructured sensory environment (as done in this work) and in leeches exposed to specific sensory stimuli or during selective intracellular stimulation of command neurons (or extracellular stimulation of the connective fibers). Another challenge will be to correlate the spontaneous activity of the identified command neurons with the overall Markov process controlling the transitions among different behaviors. The biophysical properties of these neurons and of their synapses are likely to take part in establishing the duration of different behaviors and transitions from one behavior to another. Because the duration distribution of long stationary states, swimming, and crawling episodes is on the order of tens of seconds, the dynamics of activation and deactivation of these command neurons should occur on a time scale much longer than the usual synaptic delays. Unraveling the biophysical and neurobiological mechanisms underlying such long-term neuronal modulations appears to be the key to understanding the process of behavioral decision making in the leech.

Spatial distribution during stationary states

Although the transitions among different behavioral states satisfied the Markov property of “no memory,” a significant correlation between the orientations (as well as the elongations) of the leech body in consecutive stationary states has been found (Fig. 2E-G). This may be partially accounted for by a sort of “biomechanical memory” or inertia of the leech body that, in between two close stationary states, does not have time to change its location fully; however, because the correlation has been estimated by taking into account all identified stationary states, the existence of neural circuits that effectively implement a short-term spatial memory cannot be excluded. This possibility is also supported by the fact that, in the absence of directional sensory stimuli, the distribution of leech body orientations during stationary states is isotropic (Fig. 2A-D). Therefore, the correlation between consecutive stationary orientations cannot be attributed to the tendency of the leech to rest along some preferred direction(s). On the other hand, the finding that bright flashing lights located at specific orientations can destroy the isotropy of the leech resting orientations (Fig. 3) is easily explained as the direct and continuous effect of the sensory stimuli on the decision process, without the need to hypothesize any short-term spatial memory.

Exploratory behavior

The exploratory behavior is a highly irregular, apparently unpredictable motion, however, in which we were able to identify various recurring patterns or motifs, suggestive of a deterministic dynamics (Fig. 5). Although a test for a deterministic dynamic attractor failed when applied to a long piece (1000 s) of exploratory time series, these recurrent motifs were consistently found only in the original time series but not in its phase-randomized surrogate. This suggests that nonlinear neural pattern generators could actually produce deterministic pieces of exploratory motion embedded in an overall effectively random explorative behavior. Also, the exploratory motion could simply be continuously broken by the sensory feedback provided by the external environment (which is likely to be maximal during exploration), thereby preventing the existence of a global dynamic attractor. In any case, the existence of such motifs is likely to have a deep physiological and behavioral relevance.

Despite the heterogeneity of the activated motor patterns, as suggested by the broad distribution of the head velocity (Fig. 4E,F) and by the variety of identified recurrent motifs (Fig. 5), the residence time distribution of exploratory states could be fit by a single exponential (Fig. 1D). Moreover, transition probabilities involving exploratory states satisfied the Markov first-order dependence conditions as well as stationary and swimming states (Fig. 6C). This strongly suggests that the decision process underlying initiation and termination of the exploratory behavior should share the same neural architecture that triggers or inhibits more stereotyped behaviors such as swimming or crawling. Therefore, we are tempted to speculate that a network of “exploratory” command neurons should exist in the leech head ganglion. Similar to what was found for swimming, this high-level network could simply decide when to start and stop an exploratory episode, whereas more specialized intra-ganglia and inter-ganglia circuits could work as nonlinear pattern generators of the exploratory patterns. These speculations will be the subject of future investigations.

Short pauses that adjoined exploratory states were included in the exploratory states themselves, because their statistics were clearly distinct from that of longer stationary states (Figs. 1F, 6B). Therefore, short pauses could be considered a special kind of exploratory pattern. They could represent “intermittent locomotion,” in which an animal uses brief pauses to improve the quality of sensory information to guide its behavior in response to external stimuli. This behavior has been observed in animals ranging from protozoans to mammals (Kramer and McLaughlin, 2001).

Final remarks

This work shows that the spontaneous behavior of the leech can be quantitatively analyzed and characterized for long periods of time, providing enough data for an adequate statistical analysis. Formal schemes for the behavior can then be compared with experimental data, and mathematical models can be obtained. Such an approach has already been applied successfully to the analysis of behavior of different species, from invertebrates to primates (Bressers et al., 1995; Haccou and Meelis, 1992; Chen et al., 2002); however, the leech offers the great advantage of accessibility of the nervous system to intracellular and extracellular recordings and a growing knowledge of the command interneurons responsible for shaping its behavior. Revealing the neural correlates of the Markov process and of the exploratory patterns found in the present investigation could lead to a better general quantitative theory of animal behavior relating to neuroeconomics (Glimcher, 2003) and game theory (Schultz, 2004).