INTRODUCTION

Top Abstract Introduction Methods Results Discussion References

What is the optimal way to describe the behavior of a dynamical system? This question has a special interest for neuroscientists, because the dynamics of the nervous system seems intractably complex. Although much effort has been made in recent years to characterize neuronal complexity using tools developed to decipher nonlinear systems (Basar, 1990; Skarda and Freeman, 1987; King, 1991; Garfinkel, 1983), insight resulting from such new analysis has been limited (Rapp, 1993). What we have learned is that there are deterministic dynamics present in neuronal behavior from single cells (Aihara and Matsumoto, 1986; Mpitsos et al., 1988; Hoffman et al., 1995), from ensembles of neurons (Chang et al., 1994; Schiff et al., 1994; Hayashi and Ishizuka, 1995), and even from large-scale measurements (Rapp et al., 1989; Casdagli et al., 1996; Scott and Schiff, 1995) such as the electroencephalogram (EEG), that cannot be fully described with simple linear models. We have also learned that there is synchrony between neurons that can only be detected with nonlinear measures (Schiff et al., 1996). Perhaps the most practical application of nonlinear techniques to biological systems has been dynamical control of cardiac (Garfinkel et al., 1992) and neuronal (Schiff et al., 1994b) tissues. By exploiting the natural dynamics of the system, these techniques were used to stabilize or destabilize heart beats and neuronal firing with minimum perturbation. In this paradigm, control was achieved through stabilization of an "unstable periodic orbit" (UPO) embedded within the dynamics (Ott et al., 1990).

In a mathematical space whose coordinates represent the state of a dynamical system (state space), periodic orbits are the equilibrium states. If all of the periodic orbits in this abstract dynamical landscape are unstable, the system's temporal evolution will never settle down to any one of them. Instead, the system's behavior wanders incessantly in a sequence of close approaches to these orbits. The more unstable an orbit, the less time the system spends near it. UPOs form the "skeleton" of nonlinear dynamics, and even the behavior of chaotic systems can be characterized by an infinite set of these orbits (Auerbach et al., 1987; Cvitanovic, 1988).

One can build a model of the dynamics of a system by counting and characterizing its UPOs in a hierarchy of orbits with increasing periodicity. The accuracy of such a model can be improved by progressively adding longer period orbits to the hierarchy. The dynamical landscape can then be tessellated into regions of state space centered around these UPOs (Artuso et al., 1990a,b). Orbit locations and stabilities provide short-term prediction for the future state of the system (Pawelzik and Schuster, 1991). Furthermore, if the system is nonstationary because of slow parametric variations, this can be detected through the temporal evolution of the UPOs. With a full, possibly infinite set of UPOs and their stabilities, one can calculate thermodynamic properties of a dynamical system such as entropy and dimension. This would be of little experimental relevance were it not for the finding of Cvitanovic and colleagues (Auerbach et al., 1987; Artuso et al., 1990a) that good estimates of thermodynamic properties can be obtained by using just the short orbitsthe ones most accessible experimentally. Numerous theoretical systems were shown to be well described through this approach (Artuso et al. 1990b), but until recently procedures to rigorously identify UPOs from noisy experimental data were inadequate.

In 1995, Witkowski et al. (1995) statistically confirmed the existence of unstable period-1 orbits in short biological data sets, from fibrillating dog ventricular myocardium, by comparing the detection frequency of period-1 orbits in experimental versus surrogate data, i.e., stochastic sequences generated with statistical properties similar to those of the original data (Theiler et al., 1992). A more general approach to UPO identification was offered by Pierson and Moss (1995). Their method, like Witkowski's, relied on the recurrence of patterns in state space. By using this technique, Pei and Moss (1996a) were spectacularly successful in establishing the existence of UPOs in the crayfish caudal photoreceptor. Additional research has employed this recurrence technique to identify UPOs in catfish electroreceptors (Braun et al., 1997) and teleost Mauthner cells' synaptic noise (Faure and Korn, 1997). However, there are significant shortcomings of a strict recurrence approach to UPO identification. Recurrence requires that a system's state returns repeatedly near an orbit (Lathrop and Kostelich, 1989), yet such events may be rare in finite data sets. This is further exacerbated in biological data sets, which are typically short and nonstationary. In addition, whereas the above recurrence methods principally addressed the question of existence of UPOs, they are less suited to the enumeration of distinct orbits, especially when hierarchies of orbits with higher periodicity are present.

Our group has made substantial progress in the identification of UPOs from experimental data (So et al., 1996, 1997). We developed a transformation utilizing the local dynamics of the system such that the transformed data are concentrated about distinct UPOs. The transformation acts as a dynamical lens to enhance the probability measure about the UPOs in state space. This probability enhancement helps to offset the frequent scarcity of trajectories near UPOs. In addition, we have significantly improved the ability to identify complex higher period orbits by using fragments of trajectories near those orbits (So et al., 1997). This technique overcomes the problem that, if an orbit was unstable, the system would rarely be seen to stay near the entire orbit. We also instituted statistics via surrogate data (Theiler et al., 1992) to establish confidence limits on the probability that the identified UPOs were not spurious. Our preliminary analysis demonstrated that this approach could be successfully applied to the identification of period-1 UPOs in neuronal ensemble data (So et al., 1997). In addition, a recent report using this transform technique confirmed the existence of period-1 UPOs in epileptiform activity from human cortex (Le Van Quyen et al., 1997).

We report here the first extensive application of these methods to neuronal dynamics. We examine activity from several organizational scales of neuronal structures: network behavior from small in vitro ensembles, activity of single cells within such ensembles, and large-scale activity from human cortical electroencephalographs (EEGs). We find that UPOs are present with high statistical confidence in all of these data, and that indeed they are pervasive in such neuronal activities. These findings suggest a novel means of characterizing neuronal dynamics.

	METHODS

Top Abstract Introduction Methods Results Discussion References

Data were collected from both transversely and longitudinally cut in vitro rat hippocampal slices (see Gluckman et al., 1996 for preparation details), as well as in vivo invasive EEG recordings from human epileptic patients. Single cell action potential spikes of CA1 neurons were measured by using whole-cell attached patch-clamp recordings from slices perfused with normal artificial cerebrospinal fluid (3.5 mM [K⁺]). Network burst firing activity was measured using extracellular potential recordings from the CA3 cell body layer, which was perfused with an elevated potassium level (7.5-10.5 mM [K⁺]). Digitized human EEG was collected from patients undergoing routine evaluation for epilepsy surgery that required implanted subdural or depth electrodes for medical purposes unrelated to this study. Epileptic extracellular interictal spikes were identified from the electrode closest to the epileptic focus. Because no automated method can reliably discriminate human epileptic spikes, we hand-edited these data sets for accuracy (Scott and Schiff, 1995). Institutional Review Board and Animal Research Committee approval from the Children's National Medical Center were obtained for this research. In each case, spike or burst events were identified from the recordings, and the series of interval lengths between events was used for analysis. The relationship between raw data, event times, t_n, and interevent intervals, I_n, is illustrated in Fig. 1 a.

View larger version (27K): [in this window] [in a new window]   FIGURE 1   Period-1 UPO extraction from dynamics of neuronal ensemble in a hippocampal slice (8.5 mM [K+]). (a) Extracellular field potential with event times, tn, and interburst intervals, In. (b) Sequence of interburst intervals. (c) Dynamical representation of intervals in a two-dimensional delay embedded state space (In versus In1). Colored sequences in b and c are data points that visited near the identified period-1 orbit. Gray arrows are stable and unstable directions of the orbit. (d) Density of transformed data (black line) and mean density of transformed surrogate data (green line). (e) Probability (from 0 to 1) of transformed data density being outside the distribution of maximum peaks observed from 30 transformed surrogate densities. The tallest peak, which reached above the 95% line, marks the position of a significant period-1 orbit at In = 1.04 s. (f) Sequence of points that approached along the stable and departed along the unstable direction of the period-1 orbit. One-step predicted values  of the actual experimental data , from the estimated local dynamics of the period-1 orbit. Lyapunov numbers u = 1.5 ± 0.3, s = 0.6 ± 0.2.

Our basic assumption is that there exists a significant deterministic component within the seemingly noisy activity of neurons and their ensembles, and therefore UPOs can be used to characterize the system's dynamics. The first step in our analysis is to use delay coordinate embedding (Takens, 1981; Sauer et al., 1991; Sauer, 1994) to reconstruct the underlying dynamics from our experimental data. In general, from a data sequence {x_n}, delay vectors _n in state space Z with dimension M and time lag  are given by _n = (x_n, x_n, ... , x_n(M1)). With a proper choice of the parameters M and , the geometric object defined in Z by the locus of points {_n} will provide a model that is topologically equivalent to the original dynamics that generated the data sequence.

To extract the unstable periodic orbits from reconstructed state space, a transformation based on observed local dynamics is applied to our data (So et al., 1996, 1997). The transformation concentrates the data around the UPOs. In particular, one can show that for noise-free dynamics, the probability distribution function of the transformed data will have singularities at the true periodic orbits. This reduces extraction of periodic orbits from experimental data to simply looking for peaks in the distribution of the transformed data.

The probability enhancement effect of this transformation can be illustrated by using a one-dimensional discrete-time dynamical system, x_n+1 = f(x_n), where f(x_n) is a nonlinear function that prescribes the evolution of the system state x_n. A period p orbit, x^*_p, is defined by the condition that the pth-iterated image of x^*_p is again x^*_p, i.e., where  denotes functional composition. We assume that the underlying dynamics described by the function f(x) is unknown to us, but that the local behavior of the dynamics, f'(x) = df/dx (x), can be estimated from a local least-squares fit to the experimental data. The periodic orbit transform g(x_n, ) of x_n for period-1 orbits x*:f(x*) = x* is then defined as where is a function defined by the estimated local dynamics f'(x_n) and by , an adjustable parameter of the transform. Periodic orbit transforms for higher dimension and periodicity can be analogously defined (So et al., 1997). With some algebra, it can be shown that for period-1 orbits x*, g(x*, ) = x* and dg/dx (x*, ) = 0, independent of . Therefore, the Taylor's expansion of g(x, ) near x* is approximately described by a quadratic function, g(x, )  x* + a(x  x*)², where a is a constant. After rewriting this equation, we have (g(x, )  x*)  a(x  x*)². This equation explains the probability-enhancing effect of the transform near x*. If we start with a uniformly distributed set of data points, the probability that a data point will land in a small interval with length  around x* is proportional to . Under the transformation, the same points will now be concentrated within a smaller interval with length proportional to ², thus enhancing the probability measure around x*. In practice, UPOs are found by looking for sharp peaks in the spatial distribution of the transformed data.

However, in a real experimental setting, data sets are usually contaminated by mixtures of dynamical and observational noise. In these cases, the observed sharp peaks in the distribution of transformed data are blurred into broad maxima, and can even be completely washed out by large noise. Comparison with amplitude-adjusted Fourier transform surrogates (Theiler et al., 1992) can be used to assess the dynamical significance of the observed peaks. Each surrogate is a realization of a linear stochastic model of the data formed with the same amplitude distribution and approximately the same autocorrelation function as the original data. Therefore, we do not expect UPO structure to exist in the surrogates. By using multiple realizations of the surrogate data and the statistics of extremes (Gumbel, 1958), we can estimate the probability that peaks observed from our data are statistically significant. In particular, for a specific peak observed in our transformed data density, we can estimate the probability that a maximal peak with greater amplitude could be found in the transformed surrogate densities. In our previous works (So et al., 1996, 1997) we have demonstrated the ability of our technique to extract UPOs from numerical models with both dynamical and observational noise components.

An additional problem with biological systems is their inherent nonstationarity. Such nonstationarity can be associated with parameter changes, which would be reflected in changes in the UPO structure. If the system parameters change slowly with respect to the natural time scale of the dynamics, one may expect the periodic orbit structure to be resolvable within short data windows in time (Pei and Moss, 1996b). Even so, when a system is operating near its bifurcation points, its periodic orbit structure might experience sudden changes, including orbit creation and destruction. We therefore base the following analysis on windowed data, in which the windows are chosen short enough to approximate stationary system states, but long enough to accumulate good statistics.

	RESULTS

Top Abstract Introduction Methods Results Discussion References

Statistically significant UPOs were found from ensemble, single-cell, and human recordings, as summarized in Table 1. We exhaustively searched for period-1 orbits within all of our collected data. The first two columns in Table 1 are, respectively, the total number of experiments performed and the percentage of experiments with statistically significant period-1 orbits within at least one window. In these analyses we considered observed period-1 orbits to be statistically significant only if the distribution function of the transformed data had peaks larger than 95% of maximal peaks observed from the surrogate data sets (30-100).

For the extracellular ensemble burst-firing, about half of the experiments in the lower potassium ranges (7.5-9.5 mM [K⁺]) and 90% of the experiments with higher potassium levels (10.5 mM [K⁺]) had statistically significant period-1 orbits. For the intracellular measurements from single cells within ensembles, 100% of the experiments had period-1 UPOs. These results represent the first statistical confirmation of the prevalence of UPOs in the spiking/bursting behavior of mammalian neuronal dynamics across two distinct levels of organization.

To extend our periodic orbit analysis to larger scale neuronal ensembles, we have collected intervals of interictal spikes from human EEG data. These data were collected from four different epileptic patients during the hour before the onset of a seizure. Two of the four patients' interictal spike sequences contained statistically significant period-1 UPOs. These results represent a third level of neuronal organization for which significant UPOs were observed.

In what follows we illustrate this analysis technique in detail. First, we discuss the detection of period-1 orbits from a short time window of extracellular in vitro data (Fig. 1). Next, we describe the extraction of a hierarchy of orbits, through period-3, using intracellular data (Figs. 2 and 3). We then demonstrate that the UPOs can be used to predict the observed neuronal behavior. Finally, we show how windowed analysis of long data series can be used to track nonstationary systems. This is illustrated with both intracellular (Fig. 4) and human EEG data (Fig. 5).

View larger version (60K): [in this window] [in a new window]   FIGURE 2   Hierarchy of UPOs from a single cell. (a) Sequence of interspike intervals (seconds) and short trace of raw data (zero mV corresponds to top of vertical calibration bar). (b) Data in two-dimensional delay-embedded state space. (c and d) Color-coded density plots of period-2 and period-3 transformed data revealing a family of period-2 (left graph) and period-3 orbits (right graph). Colors indicate probability (from 0 to 1) of transformed data density being outside the distribution of maximum peaks observed from 100 transformed surrogate densities. The most significant peaks, with probability greater than 95%, are shown in red. Six possible period-2 orbits, including two strongly significant ones, can be identified in c, and five possible period-3 orbits, including three strongly significant ones, can be identified in d.

View larger version (33K): [in this window] [in a new window]   FIGURE 3   Local dynamics around three representative orbits from the identified hierarchy in Fig. 2: a period-1 orbit (left column) at In = 0.069 s; a period-2 orbit (middle column) at In = 0.047 s, In = 0.076 s; and a period-3 orbit (right column) at In = 0.090 s, In = 0.026 s, In = 0.063 s. Red and green lines indicate unstable and stable directions of the UPOs, which can be identified by the intersections of the lines. (Top row) Points mapped near identified UPOs under the periodic orbit transform. The same sequences of data similarly colored are shown in Fig. 2 a. (Middle row) Deviations between pairs of representative trajectories plotted in successive time steps (connecting lines). (Bottom row) Based on estimated local dynamics of the UPOs, all sequences with good prediction (one-step error < 0.01 s) for at least two successive time steps are shown. , Initial points of these trajectories; , the subsequent iterates; +, predicted positions of the circles. Lyapunov numbers: period-1, u = 1.8 ± 0.7, s = 0.04 ± 0.04; period-2, u = 1.14 ± 0.09, s = 0.3 ± 0.2; period-3, u = 1.6 ± 0.5, u = 1.3 ± 0.2, s = 0.22 ± 0.07.

View larger version (72K): [in this window] [in a new window]   FIGURE 4   Tracking UPOs for nonstationary intracellular data. (a) Interspike interval sequence. Shaded region is the subsection of data used in Figs. 2 and 3. (b) Color-coded density plots of period-1 transformed data for 74 overlapping windows as a function of time. Colors indicate probability (from 0 to 1) of transformed data density being outside the distribution of maximum peaks from 30 transformed surrogate densities. The horizontal axis is evenly spaced between event numbers; tic marks are at 25-s intervals.

View larger version (52K): [in this window] [in a new window]   FIGURE 5   Tracking period-1 UPOs for nonstationary human EEG interspike intervals. The color map is as in Fig. 4. Three short traces (5 s) of raw EEG corresponding to different UPOs are shown. The horizontal axis is evenly spaced between event numbers; tic marks are at 200-s intervals.

Detection of UPOs

A typical example of a period-1 orbit found in the extracellular recordings is shown in Fig. 1. This figure is laid out in the same order as the steps in our analysis. First, we extract event times and interevent intervals from the raw data (Fig. 1 a). We next make a dynamical representation of the interval sequence (Fig. 1 b) in a delay embedded state space (Fig. 1 c). We then apply the period-1 transform on each data point in this state space, and compute the distribution of transformed data (Fig. 1 d, black line). For statistical comparison, surrogate data sets are also transformed, and their mean density shown (Fig. 1 d, green line). The maximal peaks from each surrogate are then extracted, and their distribution computed. Finally, the density of the transformed data is interpreted in terms of the probability of being outside the distribution of surrogate maximal peaks (Fig. 1 e). The tallest peak in Fig. 1 e, which crosses the 95% line, indicates the existence of a period-1 orbit at I_n = 1.04 s.

The colored sequences in Fig. 1, b and d, are trajectories that visited near the identified period-1 orbit. These points were chosen because, under the transformation, they map into the peak at the period-1 orbit in Fig. 1 d. Once we identify points that map into an orbit, we use those points to estimate the local dynamics near that orbit in state space. The gray arrows in Fig. 1, c and f, are estimated unstable and stable directions of the identified orbitin this case, a saddle node. In a two-dimensional state space, trajectories near an unstable saddle node will roughly approach it along the stable direction and then depart along the unstable direction. An example of this behavior is illustrated by the closed circles in Fig. 1 f.

By using the estimated local dynamics, one can predict the next point in the series based on the current one. Specifically, we use a linear fit to the local dynamics to map the current deviation from the orbit in state space to the deviation at the next iterate. Predicted forward iterates for the closed circles in Fig. 1 f are plotted as open circles. The predicted values are excellent estimates of the observed next intervals. The ability to rigorously locate UPOs and to extract their local dynamics for prediction forms the basis for control of chaotic (Ott et al., 1990) and nonchaotic (Christini and Collins, 1995) systems.

UPO hierarchy

As discussed above, UPOs form a skeletal structure for the underlying dynamics. The accuracy of this approach improves as orbits with higher periods are included. In Fig. 2, a sequence of the 512 interspike intervals (sequence, Fig. 2 a; embedding, Fig. 2 b) from an intracellular recording was used, and a family of period-1, -2, and -3 orbits was identified. The orbit locations in a two-dimensional projection of a four-dimensional state space can be identified from the color coded-density plots shown for period-2 (Fig. 2 c) and period-3 (Fig. 2 d) transformations. (In these plots, red indicates transform densities with significance greater than 95%.)

For discrete-time dynamics in a delay embedded state space, a period p orbit will comprise p individual pieces, which are related to each other through a cyclic symmetry. Therefore, in our two-dimensional space, all period-1 orbits lie along the diagonal, all period-2 orbits are pairs of points with reflection symmetry across the diagonal, and period-3 orbits are triplets of points with triangular symmetry. In Fig. 2 c, there is one strong peak along the diagonal, which indicates a significant period-1 orbit, and six pairs of points with reflection symmetry, which indicate period-2 orbits, two of which are highly significant. In Fig. 2 d, we can enumerate five possible period-3 orbits, three of which are significant. These orbits form a complex but predictable lattice of regions in which the dynamics of the neuron is approximately periodic. Although we have only enumerated orbits up to period 3, this already represents a hierarchy of 12 equilibrium states of the complex dynamics of a neuron. Similar hierarchies were found from extracelluar data at all four levels of potassium studied (7.5, 8.5, 9.5, 10.5 mM [K⁺]).

In the next section, we discuss one of each of the orbits in Fig. 2: the period-1 orbit at I_n = 0.069 s; a period-2 orbit at I_n = 0.047 s, I_n1 = 0.076 s; and a period-3 orbit at I_n = 0.090 s, I_n1 = 0.026 s, I_n2 = 0.063 s. The positions of these orbits, their stable and unstable directions (green and red lines), and the data points that mapped under transformation into their corresponding peaks in Fig. 2 are plotted in the top row of Fig. 3. The points in the top row of Fig. 3 correspond to the intervals with the same colors in Fig. 2 a.

Dynamics near UPOs

It is our goal to approximate the neuron's full dynamics using this hierarchy by partitioning the state space into regions surrounding the UPOs. The full dynamics can then be approximated by piecing together the local dynamics within these regions. For example, the two pieces of the period-2 orbit (Fig. 3, middle column) define two regions in state space where trajectories close enough to the actual period-2 orbit will, for a short time, bounce back and forth almost periodically. In terms of our biological data, the neuron appears to fire with an approximate periodicity of a two-cycle, i.e., longer intervals interspersed between shorter ones. However, because of intrinsic instability of the UPOs and the effects of noise (or new signals coming into the neuron), this apparent periodicity will end and the system's subsequent behavior will next be approximated by other orbits. The local stability of these UPOs is typically characterized by their local Lyapunov numbers (Ott, 1993). The average duration that trajectories dwell near a UPO is described by the absolute value of its largest Lyapunov number; in discrete-time dynamics, an absolute value of a Lyapunov number greater than 1 implies instability, and an absolute value less than 1 implies stability.

To approximate the full dynamics, we first need to establish that the local dynamics near these UPOs is continuous. In other words, close trajectories near a given UPO should have similar behavior. Second, we need to estimate the local dynamics and its stability near these UPOs. Last, we need to verify that the estimated local dynamics actually predict trajectories near the UPOs.

Continuity implies that interval sequences from different times (Fig. 2 a) are mapped to trajectories in state space that stay close to each other. To illustrate this point, we choose two trajectories for each UPO that start close to each other, and link their corresponding points with lines (Fig. 3, middle row). These trajectories follow each other closely on iteration.

By connecting nearby trajectories, we also graphically illustrate the dynamical stability near the UPO. For the period-1 and period-2 orbits the local dynamics can be well described by two Lyapunov numbers (Fig. 3). In both cases, the initial deviation between the two trajectories contracts along the stable direction and expands along the unstable direction. Note that for period-2 and higher, expansion and contraction must be observed after a complete cycle, e.g., iterations 1 and 3 for period-2 and 2 and 5 for period-3. In contrast to the period-1 and -2 orbits, the local dynamics around the period-3 orbit has three Lyapunov numbers with one stable and two unstable directions (two red lines, Fig. 3, right column). The period-3 plots are two-dimensional projections of a three-dimensional space. Note that differences in the number of unstable directions among the UPOs imply that the dynamics is nonhyperbolic (Dawson et al., 1994). Again, one can see that the deviations between trajectories contract and expand along the stable and unstable directions. A hierarchy of UPOs with both expanding and contracting directions is the essential ingredient for complex deterministic dynamics.

We next verify that the estimated local dynamics actually predict trajectories near the UPOs (Pawelzik and Schuster, 1991). The estimates of the dynamics near the UPO are applied to each point in the original data to predict the next point. All trajectories with good prediction (error < 0.01 s) for at least two successive time steps are plotted in the last row of Fig. 3. The first points of these trajectories are indicated by triangles, subsequent iterates by circles, and predicted positions of the circles by pluses. Prediction was excellent for points near the period-1 and period-2 orbits. The dynamical fit near the period-3 orbit was less accurate, as was the prediction. These results are the first rigorous demonstration of a dynamically meaningful hierarchy of UPOs in neuronal dynamics.

Nonstationarity: temporal evolution of UPOs

Because of nonstationary, the periodic orbit analysis was done on windowed data. Columns 3 and 4 from Table 1 are the total number of time windows partitioned from all of the data sets and the percentage of (nonoverlapping) windows in which significant period-1 orbits were found. The percentage of time windows with significant period-1 orbits was 12-28%, compared to 50-100% of experiments with significant period-1 orbits. These percentages confirm the intrinsic nonstationarity of these data. The period-1 orbits appear and disappear within the course of an experiment.

Nonstationarity is seen in the temporal evolution of period-1 orbits from the intracellular measurements shown in Fig. 4. The original interval sequence, which represents cellular discharges over 5 min of recording, is plotted in the upper panel. The hierarchy of orbits shown in Figs. 2 and 3 was extracted from the shaded region. The period-1 transformed density plot as a function of window is shown in the lower panel. The statistical significance of the transformed data density is color coded, with yellow and red indicating high significance (red, >95% confidence) for period-1 orbits. Orbits, some with high degrees of significance, were created and destroyed throughout the experiment. As in other physical systems (So et al., 1997; Carroll et al., 1992; Gluckman et al., 1997), we can characterize the nonstationarity of these systems by tracking their periodic orbit structure as a function of time.

Nonstationarity was also seen from the temporal variability of period-1 orbits from epileptic interictal spike intervals from human EEG, as shown in Fig. 5. Three different significant period-1 orbits were found during the hour before an epileptic seizure, which began just after the recording ended. Samples of raw EEG are shown from sections corresponding to each of these UPOs. As the seizure approached, the UPOs shifted to longer time intervals (0.36 s, 0.70 s, 1.23 s). Whereas the data in Fig. 5 were recorded from a neocortical frontal lobe focus, we have also identified significant UPOs from a human hippocampal temporal lobe focus (two of four patients; see Table 1).

	DISCUSSION

Top Abstract Introduction Methods Results Discussion References

These results are the first extensive application of rigorous UPO detection (So et al., 1996, 1997) to neuronal dynamics. We demonstrated that UPOs are prevalent features across several scales of neuronal organization, from in vitro single neurons to large-scale ensembles in humans. Our results, combined with evidence from others (Pei and Moss, 1996b; Le Van Quyen et al., 1997; Braun et al., 1997; Faure and Korn, 1997), suggest that complex neuronal dynamics contain significant deterministic components. In addition, these deterministic components are experimentally accessible from short biological data sets. Our findings complement the mounting experimental evidence that spike timing in neuronal systems is important (Mainen and Sejnowski, 1995; Hopfield, 1995).

Furthermore, we have demonstrated that hierarchies of UPOs can be extracted from neuronal dynamics. This is important because one can build a model that approximates the full dynamics by counting and characterizing the first few low period orbits of a hierarchy. Theoretical work (Ruelle, 1978; Artuso et al., 1990a) has established that such hierarchies can be used to estimate basic thermodynamic properties of a dynamical system, and our work suggests that this thermodynamic formalism can be applied to neuronal data.

To verify that a model based on a UPO hierarchy is valid, we demonstrated prediction of the experimental data (Pawelzik and Schuster, 1991). This type of predictive model can be used for successful parametric control of nonlinear systems, whether they are chaotic (Ott et al., 1990) or not (Christini and Collins, 1995).

A key aspect of our analysis is the ability to handle the inherent nonstationarity of biological data. When system parameters change, the skeleton of the dynamics, the UPOs, also change. Tracking parametric changes with UPOs has been accomplished in physical systems (Carroll et al., 1992; Gluckman et al., 1997); indeed, our results demonstrate this in neuronal systems at several levels of organization. Tracking is required for UPO-based control of nonstationary systems. Furthermore, tracking could be used to detect changes in system state due to intrinsic parameter variations, such as the transition to epileptic seizures, or extrinsic effects, such as electromagnetic fields Gluckman et al., 1996). In addition, recent work (Le Van Quyen et al., 1997) suggests that UPO analysis might track perceptual discrimination.

More than just forming a model for prediction and tracking, UPOs are a natural symbolic representation of a system's states. As such, we propose that UPOs form a novel symbolic language for neuronal dynamics.