- Email: [email protected]
Symbolic time-series analysis of neural data
Neurocomputing 32}33 (2000) 1073}1081
Symbolic time-series analysis of neural data夽 S. Lesher* , Li Guan, A.H. Cohen Department of Biology, Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, USA Accepted 13 January 2000
Abstract The symbolic representation of experimental data o!ers a potentially powerful tool for studying dynamic behavior and for model "tting. We show how experimental time-series data from lamprey locomotion can be mapped onto symbols. Spike time-series from bursting oscillators in the spinal cord are visualized as Poincare sections that slice through the initiation of the fast (spike) oscillations represented in a higher-dimensional phase space. This Poincare map is systematically converted to a string of symbols that capture the dynamics in the plane of slicing. The patterns in these symbol strings can be used to search for repeating behaviors, or, potentially, to "t models. 2000 Elsevier Science B.V. All rights reserved. Keywords: Timeseries; Biological oscillators; Symbolic dynamics
1. Introduction The lamprey is a segmented primitive "sh. It swims by waves of muscular contraction that move along the body. Contractions in a local group of muscles cause a region to bend towards that side. This local contraction results from networks of di!erent oscillators generating output in a group of motoneurons that emerge bundled together in the ventral root. The component neurons of the oscillators and their connections have been successfully modeled, including a single-cell multi-compartment model and a 夽
This research was supported by NIMH grant MH44809 to AHC. Present address: Department of Neuroscience, Georgetown Institute for Cognitive and Computational Sciences, Georgetown University Medical Center, 3970 Reservoir Rd. Washington, DC 20007, USA. * Corresponding author. Tel.: #1-301-405-6909; fax: #1-301-314-9358. E-mail address: [email protected] (S. Lesher). 0925-2312/00/$ - see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 5 - 2 3 1 2 ( 0 0 ) 0 0 2 8 1 - 2
1074
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
single-cell single-compartment model [26,27]. The number of types of neurons, and the total number of neurons in each half-segment are known, but the number of di!erent oscillators per segment and their spatial distribution along the cord is still undetermined. What is clear is that di!erent oscillators within a segment achieve synchrony so that all are simultaneously in spiking mode or quiet. They also maintain a "xed phase lag, of 13 per segment, between oscillators in adjacent segments along the cord [14,31]. The single-cell lamprey models successfully reproduced most aspects of the experimentally observed dynamics local to a segment. Similar models have been developed for a variety of other bursting oscillators. Such models allow detailed analyses of the dynamics of bursting and bifurcations between spiking and silent (interburst) states. One such model, for neuron R15 in Aplysia, shows how fast (spiking) and the slow (bursting) systems interact. A comparison of traditional biophysical (I}<) models with more general phase-space models clari"es the limitations of the former for representing the dynamics of Ca> and other non-voltage-dependent behaviors [2,3]. However, the majority of neurophysiological data, including that for lamprey, consists of traditional biophysical (I}<) data, and time series of spike events and amplitudes. For a single-cell burster such as Aplysia R15 (and presumably for lamprey single-cell models), it is possible to make a section through the phase space of the model, `slicinga the oscillatory trajectories for individual spike dynamics at the point where voltage reaches the threshold. This slicing generates a Poincare section through the phase space, reducing higher-dimensional trajectories to a series of lower-dimensional points corresponding to traditional time-of-spike-event data (see, for example, Fig. 1A in [1]). In the case of real lamprey motor roots, a time-of-spike-event data series represents the output of a number of di!erent motoneurons presumably re#ecting di!erent oscillators within a segment. Since their spiking is coordinated, the di!erent oscillators must be synchronized in some fashion. Spike timings re#ect the dynamics of the synchronization, as well as the dynamics of individual bursters (which in turn includes interactions between the fast (spiking) and slow (bursting) oscillators in the system). Earlier, we found substantial evidence that there is a deterministic interaction among spikes. We used the delta}epsilon method developed by Kaplan that looks for continuity in the behavior of neighboring points in phase space [13]. Under this method, in a wide variety of lamprey data, spike-to-spike events showed about as much evidence of determinism as did the much more frequently studied burst-to-burst timings [15,16]. General embedding theorems extended to interspike intervals show that, at least in theory, such timings can be used to reconstruct all of the dynamics [21,22].
Butera suggests that a symbolic dynamics based simply on counting the number of action potentials in each burst could be used to analyze a bursting time series for recurring patterns.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1075
It will take much additional work, both experimental and modeling, to be able to map interspike intervals in lamprey onto the underlying dynamics. What we show here is how spike timings can be represented abstractly as a symbol string in a way that captures some of the local dynamics, and possible future directions for this work.
2. Methods We assume that we have a trajectory in a three-dimensional phase space where one dimension is the time interval between two events, say two spikes, ¹ , a second L dimension is the next time interval, ¹ , and the third dimension lumps together all L> the Hodgkin}Huxley-like factors that change cyclically to cause spiking. Assume that this trajectory is a limit cycle with no noise, and that we look only at the twodimensional plane de"ned by ¹ and ¹ . The trajectory of the system in threeL L> dimensional phase space is reduced to a single point on the two-dimensional Poincare section. In an experimental time series of interspike (or interburst) intervals from neurophysiological systems, a stable limit cycle would have all these intervals equal, or approximately so, and a single point would appear on the diagonal at ¹ "¹ . L L> In data from most experimental systems including the lamprey, the single point of the stable limit cycles is replaced by a pattern of points as a more complex threedimensional trajectory slices through the, ¹ "¹ plane, because if ¹ is not equal L L> L to ¹ , points will not fall on the diagonal. These unstable periodic orbits produce L> characteristic sequences of points in Poincare section. But noise and non-stationarity and feedback control and inputs from other oscillators alter the trajectories in neurophysiological systems so that a given sequence of points is seldom exactly revisited. Furthermore, any speci"c pattern of points, such as the one characterizing a #ip saddle or other unstable periodic orbit, can be found in a random sequence of interspike intervals. Thus, it is necessary to apply statistical tests to see whether such sequences appear in any given data set more frequently than would be expected by chance. We scrambled the order of the interspike intervals in original time-series sequences to produce a surrogate data set, and then searched for patterns such as those typical of #ip saddles. We repeated this 20}40 times. The mean number of such patterns in this surrogate data was compared with the number in the original sequence. Using a statistic Z"(number}in}data!mean}number}in}surrogates)/std}dev}of!surrogates), we showed that there are signi"cant di!erences in the frequency of di!erent interval patterns between the original data and random rearrangements [23,28,29]. With this basic technique, unstable periodic orbits have been found in a variety of electrophysiological time series [11,18,24,32] at a rate signi"cantly higher than would be expected by chance. More generally, unstable periodic orbits may be characterized and classi"ed in terms of the properties of a local linearized system of equations. This system of equations describes the neighborhood of the unstable "xed points in the
1076
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
two-dimensional section through unstable periodic orbits [12]. A related approach has been used to analyze the local dynamics of two di!erent models of the Aplysia R15 bursting neuron at points approximately equivalent to those in our Poincare section [1}3]. For "ve intervals ¹ , we used a local linearization involving simultaneous solution G of three equations of the form ¹ "a#b¹ #c¹ , L> L> L> ¹ "a#b¹ #c¹ , L> L> L ¹ "a#b¹ #¹ L> L L\ We can solve for a, b, and c exactly. (We often used a larger window, in which case there are additional equations of this form for each additional interval. Such a system can be solved by least-squares best "t.) We generated our sequence of ¹ using a sliding window that was moved along, one G interval at a time, through the entire data series. From the solution of these equations for each successive window, this procedure generates a matrix of coe$cients, the Jacobean. The behavior of the system in this local linear region is re#ected in the Jacobean, (@A ). The values of b and c for each solution can be mapped on to a trace-determinant plane, which classi"es the behavior of such two-dimensional maps. Here b, the trace, is the ordinant, and }c, the determinant, is the abscissa. The plane is further divided by the curve of the discriminant"b!4ac, two sloping lines: trace!determinant"1, and trace#determinant"!1, and a line connecting the tangencies of the discriminant and these diagonals, at determinant"1 (see Fig. 1). Using this approach, we categorized the local behavior of the system and how it changes from window to window, where each window is `centereda on ¹ . We L did this by generating for each window a symbol representing the area of the trace}determinant plane for the [!cb] vector for the solution of equations for that window [9,10]. Fig. 1 shows how di!erent symbols: *, #, 䉭, 䉫, 䊐, *, ¤ are assigned to di!erent regions of the trace}determinant plane, representing the local behavior of systems of equations with trace and determinant mapping to these regions. Fig. 2 shows how the system dynamics, as captured in a sliding window sequence of interspike intervals, can be mapped by the solution of local linear equations relating these intervals into a sequence of symbols. This sequence captures some of the local dynamic behavior represented in the intervals. Using this method, we mapped the original time series of interspike intervals into strings of symbols. We applied the same method to categorize random sequences (surrogate data) generated from the original time series. We showed previously that the distribution of categories in the original data di!ers from that found in the surrogates [17].
Describes the use of Jacobeans from local segmentation of phase space for a somewhat di!erent purpose.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1077
Fig. 1. Trace}determinant plane representing the behavior of systems of two-dimensional maps. We make our system two-dimensional by setting ; "¹ ; iOn and substituting so that we have a pair of G> G equations ¹ "a#b¹ #c; , ; "¹ L> L L L> L for which we can "nd the Jacobean and descriptors derived from it (trace, determinant, discriminant, eigenvalues). We have arbitrarily assigned symbols: *, #, 䉭, 䉫, 䊐, *, P, to represent di!erent regions of the trace}determinant plane.
3. Results Ideally, when representing a continuous system with discrete symbols we would like to "nd an optimal coding. This is equivalent to "nding the smallest number of bits with which to represent alternative classi"cations. This optimal classi"cation is equivalent to generating a representation where all classes are equally probable. Like other optimal codes, it involves maximizing some equivalent of a Shannon-like information measure [5,25]. An optimal experimental partition can be found heuristically even when the theoretically optimal generative partition does not exist because of noise [4,6]. The coding described here is not optimal: some regions of the trace}determinant plane, especially stable foci, often appear disproportionately. We experimented with various approaches for deriving more optimal codes, but decided to stay with this Strong et al. [25] o!er an alternative approach using information theory for such data.
1078
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
Fig. 2. Diagram representing how a window of interspike intervals is slid through a time series of such intervals. Using the equations in the text (see also Fig. 1), the behavior of this window of intervals is mapped onto an area of the trace}determinant plane, and hence a symbol. The sequence of windows of interspike intervals thus becomes a sequence of symbols representing a sequence of local dynamics as captured in the interspike intervals.
representation because of its simple relation to the dynamics of point-to-point mapping in the Poincare section. Since we wanted to see quickly if there was any promise in this approach, our preliminary analysis was based on a brute force examination of various pairs and triples of our symbols, and their correlation with spiking vs. quiescent subsequences. We created symbolic data sequences, e.g. 2##*#*2 from various preparations using the methods described above. We separated the data set into two subsequences, spiking and quiescent interbursts, on the basis of histograms of the raw time intervals, using a threshold adjusted for each data set. A sequence of states, whether elements of a genome +i.e. A, T, C, G,, or a symbolic coarse-grained representation of a dynamical system such as described here, can be used to generate a set of transition probabilities transition}probability}A to B number}of}times}A is followed by B}in}selected}subsequence " number}of}times}A is followed by B}in}whole}sequence suitable for a simple Markov model for su$ciently long sequences.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1079
Table 1 Two symbols Change in interval
Transition
Symbol 1
Symbol 2
% correct
Increasing interval Decreasing interval Increasing interval Decreasing interval
SpikingPinterburst InterburstPspiking SpikingPinterburst InterburstPspiking
¤ # 䉭 #
# 䊐 # *
100 100 74 94
Table 2 Three symbols Change in interval
Transition
Symbol 1
Symbol 2
Symbol 3
% correct
Increasing interval Decreasing interval
SpikingPinterburst InterburstPspiking
䉭 #
# *
# *
85 100
We calculated the transition probabilities for all pairs of symbols, e.g. ##,#*, **, etc. separately for each group of subsequences. After eliminating data from relatively rare transitions, we looked for symbol sequences where the transition probabilities di!ered between the two subsequences. We then looked to see how often such distinguishing subsequences correlated with actual transitions between spiking and interburst intervals. We found that the two-symbol sequence `¤#a, occurred only where the time-series intervals of the two symbols increased (spikingPinterburst) The two-symbol sequence `#䊐a is always associated with decreasing time intervals (interburstPspiking). Other two-symbol sequences were partially predictive: `䉭#a was associated with spikingPinterburst 74% of the time, while `#*a was associated with interburstPspiking 93% of the time. The predictive performance could be improved by adding a third symbol: `䉭##a predicted spikingPinterburst 85% of the time, while `#**a predicted interburst P spiking 100% of the time (Tables 1 and 2). Because there is an improvement in predictive performance in going from two symbols to three, we assume that the state of the system is modeled better with the longer symbol sequence.
4. Discussion We show here how we go from a time series of interspike intervals through a linear model of local dynamics within a window of such intervals to a symbol associated with this local behavior. We show that short sequences of these symbols correlate with
1080
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
changes in state. Hence, patterns in sequences re#ect patterns (correlations) in the underlying dynamic behavior [19,20,30]. Although the symbol sequences in the tables above were very consistently associated with the changes in state (here changes in interval length), apparently similar changes in state were also associated with alternative symbol sequences. This suggests that the underlying state of the system is not completely represented by our simple set of symbols. To achieve a better mapping we will need to go to a Markov chain models with hidden states * that is, states that are not directly represented by the symbols in the output. Transition probabilities for the more complex hidden Markov models are derived using algorithms that repeatedly iterate this transition matrix, multiplying conditional probabilities. Hidden Markov models have similarities to the state-space reconstruction (the reconstruction of a trajectory from the sequence of ¹ , that is, L interevent intervals) used in time-series analysis of nonlinear dynamical systems [8]. However, we feel that it may be more productive in the immediate future to take models of lamprey bursting behavior and adjust parameters until they replicate the behavior observed in experimental spike trains. In particular, we would seek to have models replicate the observed distributions of di!erent symbol output sequences [7]. While it may be possible to achieve this with single-cell models, we think that part of the dynamics re#ected in the sequence of interspike intervals re#ects synchronization behaviors. Synchronization must occur both among oscillators producing output local to a given segment and with oscillators elsewhere along the cord. By testing the symbolic output of di!erent forms of coupling between single-cell oscillator models, and adjusting parameters until this output corresponds to the symbolic output from experimental data, it will be possible to iteratively re"ne our models of local oscillator-to-oscillator coupling.
References [1] R.J. Butera, Multirhythmic bursting, Chaos 8 (1998) 274}284. [2] R.J. Butera, J.W. Clark, J.H. Byrne, Transient responses of a modeled bursting neuron: analysis with equilibrium and averaged nullclines, Biol. Cybernet. 77 (1997) 307}322. [3] R.J. Butera, J.W. Clark, C.C. Canavier, D.A. Baxter, J.H. Byrne, Analysis of the e!ects of modulatory agents on a modeled bursting neuron: dynamic interactions between voltage and calcium dependent systems, J. Comput. Neurosci. 2 (1995) 19}44. [4] J.P. Crutch"eld, N.H. Packard, Symbolic dynamics of noisy chaos, Physica D 7 (1983) 201}223. [5] D.L. Dowe, K. Prank, Complexity and information-theoretic approaches to biology, Paci"c Symposium on Biocomputing, Vol. 3, 1998, pp. 557}558. [6] C.E.A. Finney, Bibliography of symbolic time series analysis a.k.a. symbolic-sequence analysis, 1998. http://www-chaos.utk.edu/bibSTSA.html. [7] C.E.A. Finney, J.B. Green Jr., C.S. Daw, Symbolic time-series analysis of engine combustion measurements, SAE International Congress and Exposition, Detroit, MI1998, available via [6], above. Victor Purpura [30] use an alternative method related to edit-distance formalism from genome sequence analysis (and molecular phylogeny classi"cation) to compare (measure distances between )spike trains.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1081
[8] A.M. Fraser, A. Dimitriadis, Forecasting probability densities by using hidden Markov models with mixed states, in: A.S. Weigend, N. A. (Eds.), Time Series Prediction: Forecasting the Future and Understanding the Past, Addison-Wesley, Reading, MA. [9] G. Froyland, Estimating physical invariants and space averages of dynamical systems indicators, Bull. Austral. Math. Soc 56 (1997) 157}159. [10] G. Froyland, Estimating physical invariants and space averages of dynamical systems indicators, Ph.D. Thesis, The University of Western Australia. [11] A. Gar"nkel, M.L. Spano, W.L. Ditto, J.N. Weiss, Controlling cardiac chaos, Science 257 (1992) 1230}1235. [12] M.W. Hirsch, S. Smale, Di!erential Equations, Dynamical Systems, and Linear Algebra, Academic Press, San Diego, 1974. [13] D.T. Kaplan, Exceptional events as evidence for determinism, Physica D 73 (1994) 38}48. [14] N. Kopell, G.B. Ermentrout, T.L. Williams, On chains of oscillators forced at one end. SIAM J. Appl. Math. 51 (1991) 1397}1417. [15] S. Lesher, Stable lamprey swimming on a skeleton of unstable periodic orbits, Ph.D. Dissertation, University of Maryland, 1998. [16] S. Lesher, P.E. Latham, A.H. Cohen, Evidence for deterministic dynamics in lamprey spinal cord, Soc. Neurosci. 23 (1997) 205 (Abstract 86.4).. [17] S. Lesher, M.L. Spano, N.M. Mellen, L. Guan, S. Dykstra, A.H. Cohen, Evidence for unstable periodic orbits in intact swimming lampreys isolated, spinal cords, and intermediate preparations, Ann. NY Acad. Sci. 860 (1998) 486}491. [18] D. Pierson, F. Moss, Detecting periodic unstable points in noisy chaotic and limit cycle attractors with applications to biology, Phys. Rev. Lett. 75 (1995) 2124}2127. [19] P.D. Roberts, Classi"cation of rhythmic patterns in the stomatogastric gangion, Neuroscience 81 (1997) 281}296. [20] P.D. Roberts, Classi"cation of temporal patterns in dynamic biological networks, Neural Comput. 10 (1998) 1831}1846. [21] T. Sauer, Reconstruction of dynamical systems from interspike intervals, Phys. Rev. Lett. 72 (1994) 3811. [22] T. Sauer, J.A. Yorke, M. Casdagli, Embedology, J. Stat. Phys. 65 (1991) 579}616. [23] T. Schreiber, A. Schmitz, Improved surrogate data for nonlinearity tests, Phys. Rev. Lett. 77 (1996) 635. [24] P. So, E. Ott, T. Sauer, B.J. Gluckman, C. Grebogi, S.J. Schi!, Extracting unstable periodic orbits from chaotic time series data, Phys. Rev. E 55 (1997) 5398}5417. [25] S.P. Strong, R.R. de Ruyter van Stevenick, W. Bialek, R. Koberle, On the application of information theory to neural spike trains, Paci"c Symposium on Biocomputing, Vol. 3, pp. 619-630. [26] J. TegneH r, A single-compartmental formulation and mathematical analysis of a biophysically based neuron model, Soc. Neurosci. 24 (1998) 1814 (Abstract 719.8).. [27] J. TegneH r, A. Lansner, S. Grillner, Modulation of burst frequency by calcium-dependent potassium channels in the lamprey locomotor system: dependence on the activity level, J. Comp. Neurosci. 5 (1998) 121}140. [28] J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J.D. Farmer, Testing for nonlinearity in time series: the method of surrogate data, Physica D 58 (1992) 77}94. [29] J. Theiler, D. Prichard, Using `surrogate surrogate dataa to calibrate the actual rate of false positives in tests for nonlinearity in time series, Fields Inst. Commun. 11 (1997) 99}113. [30] J.D. Victor, K.P. Purpura, Metric-space analysis of spike trains: theory, algorithms, and application, Network-Comput. Neural Systems 8 (1997) 127}164. [31] T.L. Williams, K.A. Sigvardt, N. Kopell, G.B. Ermentrout, M.P. Remler, Forcing of coupled nonlinear oscillators: studies of intersegmental coordination in the lamprey locomotor central pattern generator, J. Neurophys. 64 (1990) 862}871. [32] F. Witkowski, K.M. Kavanagh, P.A. Penkoske, R. Plonsey, M.L. Spano, W.L. Ditto, D.T. Kaplan, Evidence for determinism in ventricular "brillation, Phys. Rev. Lett. 75 (1995) 1230}1233.
Symbolic time-series analysis of neural data夽 S. Lesher* , Li Guan, A.H. Cohen Department of Biology, Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, USA Accepted 13 January 2000
Abstract The symbolic representation of experimental data o!ers a potentially powerful tool for studying dynamic behavior and for model "tting. We show how experimental time-series data from lamprey locomotion can be mapped onto symbols. Spike time-series from bursting oscillators in the spinal cord are visualized as Poincare sections that slice through the initiation of the fast (spike) oscillations represented in a higher-dimensional phase space. This Poincare map is systematically converted to a string of symbols that capture the dynamics in the plane of slicing. The patterns in these symbol strings can be used to search for repeating behaviors, or, potentially, to "t models. 2000 Elsevier Science B.V. All rights reserved. Keywords: Timeseries; Biological oscillators; Symbolic dynamics
1. Introduction The lamprey is a segmented primitive "sh. It swims by waves of muscular contraction that move along the body. Contractions in a local group of muscles cause a region to bend towards that side. This local contraction results from networks of di!erent oscillators generating output in a group of motoneurons that emerge bundled together in the ventral root. The component neurons of the oscillators and their connections have been successfully modeled, including a single-cell multi-compartment model and a 夽
This research was supported by NIMH grant MH44809 to AHC. Present address: Department of Neuroscience, Georgetown Institute for Cognitive and Computational Sciences, Georgetown University Medical Center, 3970 Reservoir Rd. Washington, DC 20007, USA. * Corresponding author. Tel.: #1-301-405-6909; fax: #1-301-314-9358. E-mail address: [email protected] (S. Lesher). 0925-2312/00/$ - see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 5 - 2 3 1 2 ( 0 0 ) 0 0 2 8 1 - 2
1074
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
single-cell single-compartment model [26,27]. The number of types of neurons, and the total number of neurons in each half-segment are known, but the number of di!erent oscillators per segment and their spatial distribution along the cord is still undetermined. What is clear is that di!erent oscillators within a segment achieve synchrony so that all are simultaneously in spiking mode or quiet. They also maintain a "xed phase lag, of 13 per segment, between oscillators in adjacent segments along the cord [14,31]. The single-cell lamprey models successfully reproduced most aspects of the experimentally observed dynamics local to a segment. Similar models have been developed for a variety of other bursting oscillators. Such models allow detailed analyses of the dynamics of bursting and bifurcations between spiking and silent (interburst) states. One such model, for neuron R15 in Aplysia, shows how fast (spiking) and the slow (bursting) systems interact. A comparison of traditional biophysical (I}<) models with more general phase-space models clari"es the limitations of the former for representing the dynamics of Ca> and other non-voltage-dependent behaviors [2,3]. However, the majority of neurophysiological data, including that for lamprey, consists of traditional biophysical (I}<) data, and time series of spike events and amplitudes. For a single-cell burster such as Aplysia R15 (and presumably for lamprey single-cell models), it is possible to make a section through the phase space of the model, `slicinga the oscillatory trajectories for individual spike dynamics at the point where voltage reaches the threshold. This slicing generates a Poincare section through the phase space, reducing higher-dimensional trajectories to a series of lower-dimensional points corresponding to traditional time-of-spike-event data (see, for example, Fig. 1A in [1]). In the case of real lamprey motor roots, a time-of-spike-event data series represents the output of a number of di!erent motoneurons presumably re#ecting di!erent oscillators within a segment. Since their spiking is coordinated, the di!erent oscillators must be synchronized in some fashion. Spike timings re#ect the dynamics of the synchronization, as well as the dynamics of individual bursters (which in turn includes interactions between the fast (spiking) and slow (bursting) oscillators in the system). Earlier, we found substantial evidence that there is a deterministic interaction among spikes. We used the delta}epsilon method developed by Kaplan that looks for continuity in the behavior of neighboring points in phase space [13]. Under this method, in a wide variety of lamprey data, spike-to-spike events showed about as much evidence of determinism as did the much more frequently studied burst-to-burst timings [15,16]. General embedding theorems extended to interspike intervals show that, at least in theory, such timings can be used to reconstruct all of the dynamics [21,22].
Butera suggests that a symbolic dynamics based simply on counting the number of action potentials in each burst could be used to analyze a bursting time series for recurring patterns.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1075
It will take much additional work, both experimental and modeling, to be able to map interspike intervals in lamprey onto the underlying dynamics. What we show here is how spike timings can be represented abstractly as a symbol string in a way that captures some of the local dynamics, and possible future directions for this work.
2. Methods We assume that we have a trajectory in a three-dimensional phase space where one dimension is the time interval between two events, say two spikes, ¹ , a second L dimension is the next time interval, ¹ , and the third dimension lumps together all L> the Hodgkin}Huxley-like factors that change cyclically to cause spiking. Assume that this trajectory is a limit cycle with no noise, and that we look only at the twodimensional plane de"ned by ¹ and ¹ . The trajectory of the system in threeL L> dimensional phase space is reduced to a single point on the two-dimensional Poincare section. In an experimental time series of interspike (or interburst) intervals from neurophysiological systems, a stable limit cycle would have all these intervals equal, or approximately so, and a single point would appear on the diagonal at ¹ "¹ . L L> In data from most experimental systems including the lamprey, the single point of the stable limit cycles is replaced by a pattern of points as a more complex threedimensional trajectory slices through the, ¹ "¹ plane, because if ¹ is not equal L L> L to ¹ , points will not fall on the diagonal. These unstable periodic orbits produce L> characteristic sequences of points in Poincare section. But noise and non-stationarity and feedback control and inputs from other oscillators alter the trajectories in neurophysiological systems so that a given sequence of points is seldom exactly revisited. Furthermore, any speci"c pattern of points, such as the one characterizing a #ip saddle or other unstable periodic orbit, can be found in a random sequence of interspike intervals. Thus, it is necessary to apply statistical tests to see whether such sequences appear in any given data set more frequently than would be expected by chance. We scrambled the order of the interspike intervals in original time-series sequences to produce a surrogate data set, and then searched for patterns such as those typical of #ip saddles. We repeated this 20}40 times. The mean number of such patterns in this surrogate data was compared with the number in the original sequence. Using a statistic Z"(number}in}data!mean}number}in}surrogates)/std}dev}of!surrogates), we showed that there are signi"cant di!erences in the frequency of di!erent interval patterns between the original data and random rearrangements [23,28,29]. With this basic technique, unstable periodic orbits have been found in a variety of electrophysiological time series [11,18,24,32] at a rate signi"cantly higher than would be expected by chance. More generally, unstable periodic orbits may be characterized and classi"ed in terms of the properties of a local linearized system of equations. This system of equations describes the neighborhood of the unstable "xed points in the
1076
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
two-dimensional section through unstable periodic orbits [12]. A related approach has been used to analyze the local dynamics of two di!erent models of the Aplysia R15 bursting neuron at points approximately equivalent to those in our Poincare section [1}3]. For "ve intervals ¹ , we used a local linearization involving simultaneous solution G of three equations of the form ¹ "a#b¹ #c¹ , L> L> L> ¹ "a#b¹ #c¹ , L> L> L ¹ "a#b¹ #¹ L> L L\ We can solve for a, b, and c exactly. (We often used a larger window, in which case there are additional equations of this form for each additional interval. Such a system can be solved by least-squares best "t.) We generated our sequence of ¹ using a sliding window that was moved along, one G interval at a time, through the entire data series. From the solution of these equations for each successive window, this procedure generates a matrix of coe$cients, the Jacobean. The behavior of the system in this local linear region is re#ected in the Jacobean, (@A ). The values of b and c for each solution can be mapped on to a trace-determinant plane, which classi"es the behavior of such two-dimensional maps. Here b, the trace, is the ordinant, and }c, the determinant, is the abscissa. The plane is further divided by the curve of the discriminant"b!4ac, two sloping lines: trace!determinant"1, and trace#determinant"!1, and a line connecting the tangencies of the discriminant and these diagonals, at determinant"1 (see Fig. 1). Using this approach, we categorized the local behavior of the system and how it changes from window to window, where each window is `centereda on ¹ . We L did this by generating for each window a symbol representing the area of the trace}determinant plane for the [!cb] vector for the solution of equations for that window [9,10]. Fig. 1 shows how di!erent symbols: *, #, 䉭, 䉫, 䊐, *, ¤ are assigned to di!erent regions of the trace}determinant plane, representing the local behavior of systems of equations with trace and determinant mapping to these regions. Fig. 2 shows how the system dynamics, as captured in a sliding window sequence of interspike intervals, can be mapped by the solution of local linear equations relating these intervals into a sequence of symbols. This sequence captures some of the local dynamic behavior represented in the intervals. Using this method, we mapped the original time series of interspike intervals into strings of symbols. We applied the same method to categorize random sequences (surrogate data) generated from the original time series. We showed previously that the distribution of categories in the original data di!ers from that found in the surrogates [17].
Describes the use of Jacobeans from local segmentation of phase space for a somewhat di!erent purpose.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1077
Fig. 1. Trace}determinant plane representing the behavior of systems of two-dimensional maps. We make our system two-dimensional by setting ; "¹ ; iOn and substituting so that we have a pair of G> G equations ¹ "a#b¹ #c; , ; "¹ L> L L L> L for which we can "nd the Jacobean and descriptors derived from it (trace, determinant, discriminant, eigenvalues). We have arbitrarily assigned symbols: *, #, 䉭, 䉫, 䊐, *, P, to represent di!erent regions of the trace}determinant plane.
3. Results Ideally, when representing a continuous system with discrete symbols we would like to "nd an optimal coding. This is equivalent to "nding the smallest number of bits with which to represent alternative classi"cations. This optimal classi"cation is equivalent to generating a representation where all classes are equally probable. Like other optimal codes, it involves maximizing some equivalent of a Shannon-like information measure [5,25]. An optimal experimental partition can be found heuristically even when the theoretically optimal generative partition does not exist because of noise [4,6]. The coding described here is not optimal: some regions of the trace}determinant plane, especially stable foci, often appear disproportionately. We experimented with various approaches for deriving more optimal codes, but decided to stay with this Strong et al. [25] o!er an alternative approach using information theory for such data.
1078
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
Fig. 2. Diagram representing how a window of interspike intervals is slid through a time series of such intervals. Using the equations in the text (see also Fig. 1), the behavior of this window of intervals is mapped onto an area of the trace}determinant plane, and hence a symbol. The sequence of windows of interspike intervals thus becomes a sequence of symbols representing a sequence of local dynamics as captured in the interspike intervals.
representation because of its simple relation to the dynamics of point-to-point mapping in the Poincare section. Since we wanted to see quickly if there was any promise in this approach, our preliminary analysis was based on a brute force examination of various pairs and triples of our symbols, and their correlation with spiking vs. quiescent subsequences. We created symbolic data sequences, e.g. 2##*#*2 from various preparations using the methods described above. We separated the data set into two subsequences, spiking and quiescent interbursts, on the basis of histograms of the raw time intervals, using a threshold adjusted for each data set. A sequence of states, whether elements of a genome +i.e. A, T, C, G,, or a symbolic coarse-grained representation of a dynamical system such as described here, can be used to generate a set of transition probabilities transition}probability}A to B number}of}times}A is followed by B}in}selected}subsequence " number}of}times}A is followed by B}in}whole}sequence suitable for a simple Markov model for su$ciently long sequences.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1079
Table 1 Two symbols Change in interval
Transition
Symbol 1
Symbol 2
% correct
Increasing interval Decreasing interval Increasing interval Decreasing interval
SpikingPinterburst InterburstPspiking SpikingPinterburst InterburstPspiking
¤ # 䉭 #
# 䊐 # *
100 100 74 94
Table 2 Three symbols Change in interval
Transition
Symbol 1
Symbol 2
Symbol 3
% correct
Increasing interval Decreasing interval
SpikingPinterburst InterburstPspiking
䉭 #
# *
# *
85 100
We calculated the transition probabilities for all pairs of symbols, e.g. ##,#*, **, etc. separately for each group of subsequences. After eliminating data from relatively rare transitions, we looked for symbol sequences where the transition probabilities di!ered between the two subsequences. We then looked to see how often such distinguishing subsequences correlated with actual transitions between spiking and interburst intervals. We found that the two-symbol sequence `¤#a, occurred only where the time-series intervals of the two symbols increased (spikingPinterburst) The two-symbol sequence `#䊐a is always associated with decreasing time intervals (interburstPspiking). Other two-symbol sequences were partially predictive: `䉭#a was associated with spikingPinterburst 74% of the time, while `#*a was associated with interburstPspiking 93% of the time. The predictive performance could be improved by adding a third symbol: `䉭##a predicted spikingPinterburst 85% of the time, while `#**a predicted interburst P spiking 100% of the time (Tables 1 and 2). Because there is an improvement in predictive performance in going from two symbols to three, we assume that the state of the system is modeled better with the longer symbol sequence.
4. Discussion We show here how we go from a time series of interspike intervals through a linear model of local dynamics within a window of such intervals to a symbol associated with this local behavior. We show that short sequences of these symbols correlate with
1080
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
changes in state. Hence, patterns in sequences re#ect patterns (correlations) in the underlying dynamic behavior [19,20,30]. Although the symbol sequences in the tables above were very consistently associated with the changes in state (here changes in interval length), apparently similar changes in state were also associated with alternative symbol sequences. This suggests that the underlying state of the system is not completely represented by our simple set of symbols. To achieve a better mapping we will need to go to a Markov chain models with hidden states * that is, states that are not directly represented by the symbols in the output. Transition probabilities for the more complex hidden Markov models are derived using algorithms that repeatedly iterate this transition matrix, multiplying conditional probabilities. Hidden Markov models have similarities to the state-space reconstruction (the reconstruction of a trajectory from the sequence of ¹ , that is, L interevent intervals) used in time-series analysis of nonlinear dynamical systems [8]. However, we feel that it may be more productive in the immediate future to take models of lamprey bursting behavior and adjust parameters until they replicate the behavior observed in experimental spike trains. In particular, we would seek to have models replicate the observed distributions of di!erent symbol output sequences [7]. While it may be possible to achieve this with single-cell models, we think that part of the dynamics re#ected in the sequence of interspike intervals re#ects synchronization behaviors. Synchronization must occur both among oscillators producing output local to a given segment and with oscillators elsewhere along the cord. By testing the symbolic output of di!erent forms of coupling between single-cell oscillator models, and adjusting parameters until this output corresponds to the symbolic output from experimental data, it will be possible to iteratively re"ne our models of local oscillator-to-oscillator coupling.
References [1] R.J. Butera, Multirhythmic bursting, Chaos 8 (1998) 274}284. [2] R.J. Butera, J.W. Clark, J.H. Byrne, Transient responses of a modeled bursting neuron: analysis with equilibrium and averaged nullclines, Biol. Cybernet. 77 (1997) 307}322. [3] R.J. Butera, J.W. Clark, C.C. Canavier, D.A. Baxter, J.H. Byrne, Analysis of the e!ects of modulatory agents on a modeled bursting neuron: dynamic interactions between voltage and calcium dependent systems, J. Comput. Neurosci. 2 (1995) 19}44. [4] J.P. Crutch"eld, N.H. Packard, Symbolic dynamics of noisy chaos, Physica D 7 (1983) 201}223. [5] D.L. Dowe, K. Prank, Complexity and information-theoretic approaches to biology, Paci"c Symposium on Biocomputing, Vol. 3, 1998, pp. 557}558. [6] C.E.A. Finney, Bibliography of symbolic time series analysis a.k.a. symbolic-sequence analysis, 1998. http://www-chaos.utk.edu/bibSTSA.html. [7] C.E.A. Finney, J.B. Green Jr., C.S. Daw, Symbolic time-series analysis of engine combustion measurements, SAE International Congress and Exposition, Detroit, MI1998, available via [6], above. Victor Purpura [30] use an alternative method related to edit-distance formalism from genome sequence analysis (and molecular phylogeny classi"cation) to compare (measure distances between )spike trains.
S. Lesher et al. / Neurocomputing 32}33 (2000) 1073}1081
1081
[8] A.M. Fraser, A. Dimitriadis, Forecasting probability densities by using hidden Markov models with mixed states, in: A.S. Weigend, N. A. (Eds.), Time Series Prediction: Forecasting the Future and Understanding the Past, Addison-Wesley, Reading, MA. [9] G. Froyland, Estimating physical invariants and space averages of dynamical systems indicators, Bull. Austral. Math. Soc 56 (1997) 157}159. [10] G. Froyland, Estimating physical invariants and space averages of dynamical systems indicators, Ph.D. Thesis, The University of Western Australia. [11] A. Gar"nkel, M.L. Spano, W.L. Ditto, J.N. Weiss, Controlling cardiac chaos, Science 257 (1992) 1230}1235. [12] M.W. Hirsch, S. Smale, Di!erential Equations, Dynamical Systems, and Linear Algebra, Academic Press, San Diego, 1974. [13] D.T. Kaplan, Exceptional events as evidence for determinism, Physica D 73 (1994) 38}48. [14] N. Kopell, G.B. Ermentrout, T.L. Williams, On chains of oscillators forced at one end. SIAM J. Appl. Math. 51 (1991) 1397}1417. [15] S. Lesher, Stable lamprey swimming on a skeleton of unstable periodic orbits, Ph.D. Dissertation, University of Maryland, 1998. [16] S. Lesher, P.E. Latham, A.H. Cohen, Evidence for deterministic dynamics in lamprey spinal cord, Soc. Neurosci. 23 (1997) 205 (Abstract 86.4).. [17] S. Lesher, M.L. Spano, N.M. Mellen, L. Guan, S. Dykstra, A.H. Cohen, Evidence for unstable periodic orbits in intact swimming lampreys isolated, spinal cords, and intermediate preparations, Ann. NY Acad. Sci. 860 (1998) 486}491. [18] D. Pierson, F. Moss, Detecting periodic unstable points in noisy chaotic and limit cycle attractors with applications to biology, Phys. Rev. Lett. 75 (1995) 2124}2127. [19] P.D. Roberts, Classi"cation of rhythmic patterns in the stomatogastric gangion, Neuroscience 81 (1997) 281}296. [20] P.D. Roberts, Classi"cation of temporal patterns in dynamic biological networks, Neural Comput. 10 (1998) 1831}1846. [21] T. Sauer, Reconstruction of dynamical systems from interspike intervals, Phys. Rev. Lett. 72 (1994) 3811. [22] T. Sauer, J.A. Yorke, M. Casdagli, Embedology, J. Stat. Phys. 65 (1991) 579}616. [23] T. Schreiber, A. Schmitz, Improved surrogate data for nonlinearity tests, Phys. Rev. Lett. 77 (1996) 635. [24] P. So, E. Ott, T. Sauer, B.J. Gluckman, C. Grebogi, S.J. Schi!, Extracting unstable periodic orbits from chaotic time series data, Phys. Rev. E 55 (1997) 5398}5417. [25] S.P. Strong, R.R. de Ruyter van Stevenick, W. Bialek, R. Koberle, On the application of information theory to neural spike trains, Paci"c Symposium on Biocomputing, Vol. 3, pp. 619-630. [26] J. TegneH r, A single-compartmental formulation and mathematical analysis of a biophysically based neuron model, Soc. Neurosci. 24 (1998) 1814 (Abstract 719.8).. [27] J. TegneH r, A. Lansner, S. Grillner, Modulation of burst frequency by calcium-dependent potassium channels in the lamprey locomotor system: dependence on the activity level, J. Comp. Neurosci. 5 (1998) 121}140. [28] J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J.D. Farmer, Testing for nonlinearity in time series: the method of surrogate data, Physica D 58 (1992) 77}94. [29] J. Theiler, D. Prichard, Using `surrogate surrogate dataa to calibrate the actual rate of false positives in tests for nonlinearity in time series, Fields Inst. Commun. 11 (1997) 99}113. [30] J.D. Victor, K.P. Purpura, Metric-space analysis of spike trains: theory, algorithms, and application, Network-Comput. Neural Systems 8 (1997) 127}164. [31] T.L. Williams, K.A. Sigvardt, N. Kopell, G.B. Ermentrout, M.P. Remler, Forcing of coupled nonlinear oscillators: studies of intersegmental coordination in the lamprey locomotor central pattern generator, J. Neurophys. 64 (1990) 862}871. [32] F. Witkowski, K.M. Kavanagh, P.A. Penkoske, R. Plonsey, M.L. Spano, W.L. Ditto, D.T. Kaplan, Evidence for determinism in ventricular "brillation, Phys. Rev. Lett. 75 (1995) 1230}1233.