Brain-Machine Interactions for Assessing the Dynamics of Neural Systems

A critical advance for brain–machine interfaces is the establishment of bi-directional communications between the nervous system and external devices. However, the signals generated by a population of neurons are expected to depend in a complex way upon poorly understood neural dynamics. We report a new technique for the identification of the dynamics of a neural population engaged in a bi-directional interaction with an external device. We placed in vitro preparations from the lamprey brainstem in a closed-loop interaction with simulated dynamical devices having different numbers of degrees of freedom. We used the observed behaviors of this composite system to assess how many independent parameters − or state variables − determine at each instant the output of the neural system. This information, known as the dynamical dimension of a system, allows predicting future behaviors based on the present state and the future inputs. A relevant novelty in this approach is the possibility to assess a computational property – the dynamical dimension of a neuronal population – through a simple experimental technique based on the bi-directional interaction with simulated dynamical devices. We present a set of results that demonstrate the possibility of obtaining stable and reliable measures of the dynamical dimension of a neural preparation.


INTRODUCTION
We report a new approach to identifying the dynamical properties of neural tissue engaged in a bi-directional interaction with an external device. We developed a brain-machine interface (BMI) system, where the dissected brain of lamprey larvae was connected bi-directionally with an external dynamical system (EDS). While BMIs are mostly applied to clinical conditions, such as complete paralysis or loss of some senses, here we consider an application to the study of the computational properties of neural systems. One can distinguish between (i) motor BMIs, where a direct communication pathway is established from the nervous system to an external device (Donoghue, 2008;Hochberg et al., 2006); (ii) sensory BMIs, where a direct communication pathway is established from an external device, e.g. an artifi cial retina, to the nervous system (Counter, 2008;Dowling, 2005;Mokwa, 2007); and (iii) bidirectional BMIs, where a bi-directional direct communication is established between the nervous system and an external device (Bakkum et al., 2007;Karniel et al., 2005;Martinoia et al., 2004). The setup describe in this paper is an example of a bi-directional BMI.
The sensory information about the current state of the external system was conveyed to the neural preparation by means of electrical pulse trains of variable frequencies, delivered by tungsten electrode. Population spike trains extracted in real time from extracellular multiunit recording were transformed into a control signal used to drive the external system. In engineering terms, this the neural preparation. The closed-loop interaction with the device facilitates the identifi cation of specifi c computational features, such as the dynamical dimension, through the creation of a BMI system capable of generating autonomous behaviors. These are behaviors that do not require the application of external inputs, whose design involves a signifi cant amount of complexity and arbitrariness in more traditional system identifi cation methods.

GENERAL PARADIGM
We consider the bi-directional signal exchange between a neural system and an external device as a means to get at the neural preparation's dynamical dimension. We measured the dynamical dimension of the preparation embedded within the closed-loop BMI system shown schematically in Figure 1. The term "neural preparation" refers to the part of the nervous system that affects the recorded signal. In general mathematical terms, the external system is described by the state and output equations: Here, x is the state of the external device, u is a control signal, y is the read-out signal, h and p are arbitrary functions, and the subscripts refer to discrete sampling time.
One can describe in the same way the neural preparation as a dynamical system. The main assumption of our study is that there exists a state representation (s) of the neural preparation such that the changes of state are completely determined by the state itself and by the input to the neural preparation: Here, i and o refer to the input and output of the neural preparation. Both s and x are vectors. In principle, there may be infi nitely many equally valid representations for s. However, its dimension (the dynamical dimension), that is the number of independent components in any numerical representation of s, is a well-defi ned property of the preparation, which we seek to determine.
The neural preparation and the external device are interconnected in a closed-loop through two instantaneous mappings: By combining these relations with the dynamics of both the external system and the neural preparation (Eqs 1-3) one obtains a single equation: where q = [x, s] T is the state of the composite system and is a composition of the states of the two sub-systems, the external device and the neural preparation, and: This equation does not hold true in the general case. It can be violated due to degeneration of the composite system into lower order. The assessment of the dynamical dimension of the neural component can also be adversely affected by (i) signifi cant difference between natural frequencies of the neural component and the external system, and (ii) by input signals driving the neural component into saturation. We have analyzed each of these possible causes and investigate how to avoid each of them. An extensive discussion of this issue can be found in the Supplementary Material.
By construction, the behavior of the composite system does not depend on any external input, but only on its current state: this is the defi ning property of an "autonomous system." In any practical case, autonomy must be understood as an approximation, because uncontrolled external inputs cannot be entirely discounted. We assessed the dynamical dimension of the composite system, dim(q), by collecting multiple trajectories of the external device. Then, we computed the dynamical dimension of the neural preparation, dim(s), from Eq. 6, by subtracting the known dimension of the external system, dim(x), from the estimated dimension of the combined system: The main element of novelty with this approach consists in the possibility of using external systems with different dimensions. Then, it is possible to use the known difference between the dimensions of the external systems as condition to validate the estimate of the dimension of the neural preparation. This provides us with a critical tool, as the estimate of a system's dimension is typically affected by noise and by arbitrariness in some data-processing parameters.

THE NEURAL PREPARATION AND THE EXPERIMENTAL SETUP
The whole brain of the Sea Lamprey was dissected and maintained in continuously superfused Ringer's solution (see details in Alford  al., 1995;Karniel et al., 2005). The dissected neural tissue was maintained alive for up to several days. Figure 2 shows the preparation in the recording chamber. The stimulating electrode was placed among the axons of the vestibulo-reticular pathway, which originates from one of the posterior octavomotor nuclei (nOMP). This nucleus receives inputs from the vestibular capsule and its axons form synapses with the contralateral rhombencephalic neurons.
The recording electrode was placed extracellularly among the axons of the (contralateral) reticulo-spinal pathway, originating from the posterior rhombencephalic reticular nucleus (PRRN). This nucleus received monosynaptic projections from the contralateral nOMP. In behaving animals, the PRRN conveys sensory information of different nature (visual, vestibular, tactile) and central commands to the motor centers of the spinal cord (Rovainen, 1979a,b). This placement of the stimulation and recording electrodes induces predominantly excitatory responses (Alford et al., 1995).
The recorded signals were acquired at 10 kHz by a data acquisition board (National Instruments PCI-MIO-16E4) on a Pentium III 2GHz computer (Dell Computer Corporation). The recording/ stimulating paths from the workstation to the tissue followed the subsequent scheme. The recording electrode (a glass micropipette fi lled with 1 M NaCl, impedance range 1-10 MΩ) was connected to a pre-amplifi er within a Faraday cage which was in turn connected to an external amplifi er. The stimulator (Programmable Master 8 system by A.M.P.I.) delivered the stimulus pattern to the neural tissue through an isolator unit (Iso-Flex). In our experiments the stimulation was delivered through a tungsten microelectrode (impedance range of 1-2 MΩ) in the form of monophasic pulses of 200 µs duration and magnitude of 60 µA. The timing and frequency of the pulses depended on the information coming from the external device (see explanation below). Figure 3 shows a sample segment of the raw signal acquired by the recording electrode at the rate of 10 kHz during the closed-loop experiments.

SPIKE DETECTION AND ARTIFACT CANCELLATION
The procedure used for online spike detection during the closedloop experiments is described in Table 1. The minimal spike magnitude and the maximal spike duration values used for different experiments session are shown in Table 2. To avoid confusion between the stimulation artifacts and spikes, the acquired raw signal immediately following each stimulation pulse was discarded. The duration of discarded signal (the artifact cancellation period) was set to 3 ms for all experiments. The spike detection procedure in Table 1 skipped segments of the raw signal during the artifact cancellation periods.

THE EXTERNAL DYNAMICAL SYSTEMS
We used two different external systems with dynamical dimensions of two and four referred to as the 2D and the 4D system. These systems are real-time simulations of, respectively, a single point-mass and two masses connected by a spring. Both masses move along a line. For both systems, the control signal u determines the external force acting on mass m 1 , and the position of mass m 1 determines the read-out signal y. The point-masses move within a potential force  chosen because it is a monotonic and mildly non-linear function. The same expression was used in earlier studies where the neural tissue was connected to a small mobile robot (Fleming et al., 2000). The output interface was only used for normalization purposes. Different neural preparations provide different ranges of population spiking rate. The output interface normalized the spiking rate, so that the maximal observed value was mapped into 1; and converted the output signal, o, into a control signal for the external component, u, (see Figure 1): The bias, B, and the gain, C, are set to −0.2 and 10 respectively. The maximal observed output value, o max , was estimated for each experimental session in a set of preliminary runs, which are not used to collect data. Particular values of o max used in the experiments are shown in Table 2.

THE EXPERIMENTAL PROTOCOL
Each experiment consisted of a sequence of episodes, in which the BMI system operated for 20 s.
The two external systems were used in alternating episodes. A total of 20 episodes were run in each session. For each episode, the initial confi guration of the corresponding external system (positions and velocities of the point-masses) was drawn at random. The episodes were interleaved by periods of rest of 40-100 s. We found empirically that periods of this duration were suffi cient to insure stable neural responses across the experiment. We assume that different preparations might require different resting periods.
To reconstruct the dynamical dimension of the composite system we used the trajectories of the external device and not the neural activity, because the latter, being a sequence of spike counts in short-time intervals, is quite noisy, bearing the entire dimension assessment impossible. In the described setup, the external device functions as a fi lter smoothing out the noisy neural behavior. Since Table 1 loop on a segment of the recorded signal, r(t i ), starting from t 0 fi nd the next local maximum, if r(t i ) − r(t j ) > min spike magnitude and t j − t i < max spike duration Assign spike at time t i end loop fi eld, and there is a damping element between each point-mass and the ground. These dynamical features were established to induce the combined system to explore a suffi ciently wide variety of behaviors.
The Supplementary Material provides a detailed description of the two simulated dynamical systems.

THE INTERFACES
The neural tissue was connected to an EDS, thus forming a closedloop composite system. The composite system had two signal pathways: (i) a stimulating electrode conveyed the read-out signal of the EDS to the neural tissue; (ii) the signal recorded from the neural tissue determined a control signal for the EDS. The read-out signal, y, of the external system (refer to Eq. 1) was: The input interface implemented a non-linear mapping from the external system's read-out signal, y, into the input signal, i (Figure 4): The position of the simulated point-mass, x 1 , and hence the readout signal, y = x 1 , was constrained to stay within the interval [−1, 1]. Accordingly, the input signal, i, remained within [0, 1]. The "shape" parameter, a, regulated how "curved" was the function (Eq. 9); we used a = 5.This particular form of the input interface mapping was we are estimating the dynamical dimension of the entire system, it is reasonable to analyze the external device trajectories instead of applying another, additional fi lter to the neural activity. We would like to mention in this regard that the musculo-skeletal system in general (not only of Sea Lamprey) functions as a fi lter that smoothes out noisy neural activity making the entire system feasible and easily controllable. This could be a factor that determines the fi ltering properties of muscle tissue. The temporal structure of the closed-loop system implementation is shown schematically in Figure 5. A data acquisition card sampled continuously the signal from the recording electrode at 10 kHz rate. The rest of the processing took place in a periodic manner; the experiment loop was repeated once every fi xed time period, Δt = 50 ms. The CPU of the computer performed the following operations, without interrupting the data acquisition card: 1. Detect spikes from the 10 kHz signal. 2. Estimate spike rate, o by counting the number of spikes in the interval Δt. 3. Compute the output interface mapping, u(o) (Eq. 10). 4. Simulate one time-step of the EDS, i.e. compute (x 1 , x 2 , 24 of the Supplementary Material). 5. Compute the input interface mapping, i(y) (Eq. 9). 6. Draw the binary stimulation signal (BSS) using the input magnitude, i. Drawing the BSS means choosing BSS equal to 0 or 1 according to the target probability defi ned as follows. The BSS is a binary signal that determines if the pulse is emitted at the current sweep of the loop. When BSS is equal to 1, the pulse is emitted; otherwise, the pulse is not emitted. The probability of emitting a pulse is therefore: FIGURE 4 | Input interface mapping (Eq. 9). This particular shape of the input interface mapping was chosen to represent the light intensity signal as a function of robot orientation in the experiments that use the same brain-machine interface setup with a mobile robot as the external device (Fleming et al., 2000). The input signal, i, determined the current target frequency up to the maximal stimulation frequency, f max , achieved when i = 1. Particular values used in the experiments for f max are shown in Table 2. 7. If the drawn BSS is equal to 1, a command is sent to the stimulation hardware to emit one pulse.
When the elapsed time of the current sweep exceeds Δt, a new sweep is started, a new portion of data is loaded from the data acquisition card, and all the steps described above are repeated again.

COMPUTING THE DYNAMICAL DIMENSION
The state trajectories of an autonomous system do not intersect themselves or each other (Arnold, 1978). In contrast, the trajectories observed in our experiments exhibit multiple self-intersections. This may be attributed to the fact that they are one-dimensional projections from a higher dimensional state space. To reconstruct the dimension of the composite system, we unfold the observed trajectories into spaces of increasing dimensions and check for the presence of intersections (Kaplan, 1994). The lowest dimension where the trajectories do not intersect is judged to be the adequate dynamical dimension of the composite system (Abarbanel, 1996;Kaplan, 1993). To unfold trajectories into a space of a given dimension d, we used the delayed coordinates space embedding (Mane, 1981;Takens, 1981) where k enumerates the state-points along the embedded state trajectory and τ is the integer lag parameter.
To select a proper lag τ, we employed an approach that uses the fi rst (corresponding to the smallest τ) local minimum of the mutual information between samples (Abarbanel, 1996), as described in the Supplementary Material. An alternative approach for selecting the lag is to set it equal to the fi rst zero crossing of the autocorrelation of the analyzed signal. However, this only makes sense if the underlying dynamics is linear. The mutual information criterion does not have this limitation, so we used it for the data analysis.
Checking for intersections of the trajectories is a complicated task, which has engendered some controversy (Ruelle, 1990). Since the sampled trajectories are sets of measure 0 within the embedding space, one can -in principle -never observe an actual intersection (of course for a digitized signal with limited numerical resolution a probability of observing two overlapping samples is non-zero but negligibly small). In spaces of dimension >2, it is only possible to determine that two trajectory segments are close to each other. Most of the existing methods for dimension reconstruction address this diffi culty (Kaplan and Glass, 1992;Kennel et al., 1992). We propose to take advantage of a priori knowledge about the external device connected to the neural element to remove some of the uncertainties. We used the "delta-epsilon method" developed by Kaplan (Kaplan, 1994). For each dimension d, the trajectories were unfolded in d-dimensional space and pairs of consecutive points along the unfolded trajectories were analyzed. For two pairs (v i , v i + 1 ) and (v j , v j + 1 ), the distance between the fi rst points of each pair is denoted by δ: and the distance between the second points of each pair is denoted by ε: Delta-epsilon combinations with small deltas (approaching the limit δ→0) and large epsilons are interpreted as intersections of the trajectories due to projection from a higher dimensional state space onto a space of insuffi cient dimension. If small deltas always bear small epsilons, it is interpreted as an indication that δ→0 implies ε→0 that corresponds to the elimination of all intersections of the trajectories. In this case, the corresponding dimension is deemed to be an adequate estimate for the dimension of the state space.
We unfolded the observed trajectories sequentially into progressively higher dimensions using delayed coordinates space embedding (Eq. 12). For each dimension, we computed the delta-epsilon combinations and considered n combinations with the smallest deltas. To assess robustness of the method, we have varied n. To test for trajectory intersections in a given dimension d, we compared the magnitudes: for different d's (see upper panels in Figure 8 in Section "Results").
Here, Z n is the set of n delta-epsilon combinations with smallest deltas. As the dimension d increased, the corresponding magnitudes ε d decreased almost monotonically until they reached a plateau. This saturation points to the dynamical dimension of the composite system. To determine this transition in different preparations, we normalized ε d : The normalization was set so that for each system ε d spanned magnitudes between 0 and 1 (see lower panels in Figure 8 in the Section "Results"). The assessed dynamical dimension d* was established by setting a threshold h and determining the fi rst threshold crossing in the function ε d vs. d:

SURROGATE DATA ANALYSIS
To address the noise issue in the dimension analysis (see Abarbanel, 1996 for a general discussion), surrogate trajectories were generated using randomized phases. They do not refl ect any implicit dependencies and are used as examples of "pure noise." Any dimension analysis applied to these trajectories should produce the estimated dimension equal to infi nity. To generate the surrogate trajectories, a Fourier image was computed for each original trajectory. This was done after sub-sampling the trajectory with the time lag τ derived from the average mutual information measure as described above. The phases of the Fourier components were then replaced by random numbers drawn uniformly from [−π, π], while the absolute magnitudes of the Fourier

ROBUSTNESS ANALYSIS
As it is common in this type of procedure, the outcome depends on arbitrary parameters, such as the threshold, h, and the number of delta-epsilon combinations, n (Eqs 15 and 17). More defi nite and robust dimension estimates can be obtained by using two alternative external systems interconnected with the same neural preparation. We tested multiple values for h and n. For each preparation, we computed the dimension of the composite system. The corresponding differences between the estimates for the two external systems were computed for all parameter combinations (see Figure 9 in the Section "Results"). The known difference between the external systems' dimensions (two) was used to eliminate the parameter combinations that produced inconsistent results. Then, the parameter combinations were s elected that produced the most robust estimates, i.e. the largest number of equal entries in the tables.

RESULTS
To assess the dynamical dimension of the composite system, we used the read-out signal y t of the external system. For both the 2D and the 4D external systems, the read-out signal is the position of the fi rst point-mass, y = x 1 (see Eq. 8). Trajectories y t collected during the experiments and used for the subsequent analysis are shown in Figure 6. For each experimentation session, 10 trajectories were obtained with each of the two EDS (the 2D system and the 4D system). The duration of each trajectory was 20 s. The dynamical dimension was assessed on the trajectories sampled with lag τ computed individually for each set of trajectories using the mutual information measure described in the Supplementary Material. Each set of trajectories consists of all trajectories collected for one preparation using one particular external system. The average mutual information as a function of the lag τ for the trajectories obtained using three different preparations combined with two alternative EDS is shown in Figure 7. The smallest τ where the mutual information function reached a local minimum was used as the time lag for the dimension analysis. The corresponding lag is marked as τ 0 in each panel in Figure 7. The lags determined using this method lie between 8 and 10. Lags computed using the fi rst zero crossing of the autocorrelation of the analyzed signal lie between 2 and 4 for most trajectories.
To assess the dynamical dimension of the composite system, we have computed all delta-epsilon pairs (see Eqs 13 and 14) for each set of trajectories for each dimension from 1 to 20.
To compute the dynamical dimension from the delta-epsilon pairs, we have analyzed ε d (Eq. 15) and ε d (Eq. 16) as functions of dimension d. These dependencies are shown in Figure 8. The two upper panels in Figures 8A-C show ε d (Eq. 15) for different d's for the two external systems.
To address the noise issue in the dimension analysis, we have computed a surrogate trajectory for each original one. The results of the delta-epsilon analysis on the surrogate trajectories are shown in Figure 8 by red plots in upper panels in A, B, and C. These plots stand in perfect agreement with the assumption that noise has infi nite dynamical dimension. They serve as a good baseline to compare the ε d plots of the original trajectories to judge their decay and saturation.
The two lower panels in Figures 8A-C show the normalized version of ε d , ε d (Eq. 16), for different d's for the two external systems. For each preparation and for any EDS, ε d spans the interval between 0 and 1. The dotted line in each lower panel in Figures 8A-C shows the threshold level, h, used to compute the dynamical dimension d* (Eq. 17). The latter is marked in each of the lower panels in Figures 8A-C as well.
As mentioned above in Section "Materials and Methods," we must address the arbitrariness in the parameter choice, which might  Figure 9. The dynamical dimension for each preparation was assessed using two factors. First, the consistency between the results for the 2D and the 4D external systems (Figures 9A-C, left and the middle panels), which requires that the difference between the corresponding entries should be equal to 2. And, second, the general robustness with respect to parameters' values ( Figures 9A-C, right panels), i.e. the largest number of equal entries in the tables (refer to Section "Robustness Analysis"). The results of this procedure are reported in the last line of Table 2. The robustness analysis shows that (i) the dynamical dimension assessment is only mildly affected by parameter selection, and that (ii) there are large areas of parameters combinations where the assessments using the two different external systems are mutually consistent.  Eqs 15 and 16). For each preparation, the left panels show the data collected using the 2D external dynamical system, the right panels show the data collected using the 4D external dynamical system. The red plot in the upper panels show ε d computed for surrogate trajectories obtained by randomizing phases of the original trajectories. The dotted line in the lower panels shows the threshold h used to assess the dynamical dimension d*.
To investigate stability of the assessed dynamical dimension over time, we have conducted two additional experimental sessions with one of the preparations (preparation 3). These additional sessions took place after a signifi cant time delay following the fi rst session (see details in Table 2). Figure 10 shows ε d and ε d (Eqs 15 and 16) as functions of d for all three sessions conducted using preparation 3. Figure 11 shows the results of the robustness analysis with respect to the parameters n and h (Eqs 15 and 17) for all three sessions conducted using preparation 3. Figures 10 and  11 use the same representations as Figures 8 and 9 respectively. Note that Figure 10A shows the same data as Figure 8C: data from preparation 3 session 1. The same is true regarding Figures   11A and 9C. The results of all these assessments are summarized in Table 2. The dynamical dimensions estimated from the three sessions are the same. This type of analysis is needed to insure (i) that the dynamical dimension is a fi xed and well-defi ned property of a particular preparation, and (ii) that the dynamical dimension is not altered by the process used to assess it. These statements refer to long-time changes that take place (or rather do not take place) in the order of magnitude of hours; this is typical duration of one experimental session for the dynamical dimension assessment. Short-time changes caused by the probing itself cannot be In each panel, the columns correspond to different magnitudes of the threshold h (Eq. 17), and the rows correspond to different magnitudes of n (Eq. 15). The three panels in the left column show the assessed dimension, d*, for the data collected using the 2D external system. The three panels in the middle column show the assessed dimension for the data collected using the 4D external system. The three panels in the right column show the difference between the two assessments for each parameter combination. The three panels in each row show the results for a particular preparation, as marked on the left.

DISCUSSION
This paper describes a methodology for applying dynamical systems theory to the study of properties of neural tissue. We have demonstrated the feasibility of a method for estimating the dynamical dimension of a neural component in a composite neuro-robotic system. We have tested this approach in vitro, using lamprey's brain stem interconnected with a dynamical system simulated by a computer. These experiments reveal critical issues and bottlenecks of the proposed methods, as well as their weak points and limitations. The general paradigm used in this work does not try to reveal the dynamics of the neuronal network that exist independently of the external system. This work presents an approach of how to estimate the neuronal network dynamics in the context of closedloop interactions.
Our fi ndings depends on several contingent factors, including the neural population under study, the electrode placement and the selected parameter settings, such as magnitude, duration, and frequency range of the stimulation pulses. The neural preparation might involve only one neuronal population projecting monosynaptically onto another. Or, it might include multiple synaptic interactions, taking place between stimulation and recording electrodes. Furthermore, the dynamical properties of the neural preparation may be accounted for by changing concentrations of particular ions or substances at different locations along the axons that were stimulated or recorded from, or by various cellular mechanisms, such as channel and autoreceptor dynamics (Schwartz and Alford, 2000). One would expect to observe different dynamical properties under these different situations.
A specifi c conclusion regarding the experimental examples provided here is that the estimated neural dynamics may depend quite critically on the placement of the stimulating and recording electrodes. Small variations in electrode position caused a change in the observed dynamical dimension. This may happen due to any of the above-mentioned factors that affect the observed dynamics. For example, the dynamical dimension of preparations 1 and 2 was estimated to be 5, while the dynamical dimension of preparation 3 was estimated to be 3. The dependence of population signals upon the placement of the recording electrodes is a general challenge for BMI systems, and this one is no exception. This is to a large extent an unavoidable problem because the electrode placement -in our case both for stimulation and for recording -defi nes the biological network under consideration. Nevertheless, once a stable electrode placement has been attained in a BMI system, we can envision the use of this or similar techniques for determining the dynamical properties of the neural tissue that is assigned the control of an external device. While at this time this study has a limited value, we expect that the probing of neural behavior through the interaction with dynamical modules with specifi c properties -such as damping, inertia, bandwidth, etc.-will provide neural engineers with a simple way for fi ne tuning the properties of the interface based on solid and easily implementable engineering principles.
The methodology presented in this paper provides a means for assessing the dynamical dimension based only on signal properties. For its most effective use, it is recommended that it be used in combination with other physiological and anatomical techniques providing insight about the cellular structure and mechanisms that determine the observed dynamics. Here, we refer, for example, to the possibility to probe the dynamical properties of portions of neural tissue -either in vivo or in vitro -in conjunction with the application of pharmacological agents that neural plasticity mediated by nMDA or metabotropic glutamate receptors. This would provide focused information on the effects of pharmacological manipulations on the dynamical properties of biological neural networks. The fi eld of BMIs has been developing very actively for the last few years (Lebedev and Nicolelis, 2006;Martinoia et al., 2004;Novellino et al., 2007;Serruya et al., 2002;Taylor et al., 2002;Wessberg et al., 2000). The direct interaction between the nervous system and external devices has a broad range of applications. In particular, there is a promising potential to use BMIs for reconnecting with the external world patients suffering from a variety of disabling neurological disorders (Donoghue, 2002;Hochberg et al., 2006;Kennedy and Bakay, 1998;Nicolelis, 2003). The novel feature of our method is the use of the interface with an external device as a research tool. The external device is used to reveal certain characteristics of the neural preparation. In particular, we used two different external devices alternatively with the same preparation in order to address the arbitrariness in the parameter choice for data analysis. At the same time, the external device affects the behavior of the entire BMI system and implicitly shapes the neural activity of the preparation. Evidently, activity patterns exhibited by the preparation in closed-loop interaction with a particular artifi cial device may differ from those exhibited in natural conditions. In the described experiments, we do not try to simulate the natural behavior or natural neural dynamics. Therefore, all the measurements and conclusions regarding the neural preparation are strictly valid only within the context of a particular BMI system.
We need to stress that in this initial test of our procedure we have used external devices that fi ltered out high frequency components of the neural signal, so that the resulting behavior was relatively smooth. This however, needs not to be a constraint in other implementations of this approach. We have decided here to introduce an amount of low-pass fi ltering, in analogy with the smoothing properties of the musculo-skeletal system together with the spinal circuitry (Partridge, 1966(Partridge, , 1973. In the general case, different external systems excite different components of the neural dynamics and thus can be used as investigation tools aimed at different operational conditions. Ultimately, the intrinsic dimension of the neural system is likely to be higher than the estimate by our method. This is not necessarily a disadvantage, as the coupled device can deliberately be designed to fi lter out frequency components and dynamical behaviors that may be detrimental to the effective control of external devices, such as wheelchairs or robotic arms.
Our results indicate that the neural dynamics are captured by relatively simple descriptive models. Here, one must take into account that we used an extracellular electrode. This is not an implicit constraint of the method, as the same concepts can in principle be applied in combination with intracellular recording. As the task of separating single units from extracellular recordings is based on some arbitrary assumptions, we did not make any attempt to this end. We conclude with three considerations on assessing the dynamical dimension of neural preparations in experimental and theoretical investigations: 1. To design effective computational models of a neural preparation, one must take into consideration its dynamical dimension.
The dynamical dimension provides the lower bound for the number of tunable parameters for any adequate model. If the state of the neural preparation has four independent components, as we have assessed using the model-less method presented here, then any model trying to match experimental data for a similar input/output arrangement must use at least four independent parameters. 2. Our fi ndings on dynamical dimension suggest the presence of structural features of the preparation, such as connectivity patterns and dynamical properties of particular channels. These fi ndings suggest further investigations with other means. For example, pharmacological agents are used for altering channel properties (Grillner, 2003), for blocking selected ion channels (Cangiano et al., 2002), or for altering properties of multisynaptic transmission (Schwartz et al., 2005;Takahashi et al., 2001). The dynamical dimension assessment can be used in combination with these agents to provide a direct evaluation of their effect on neural dynamics.
The general framework of the present study emphasizes the interdisciplinary approach involving experimental neurobiology, engineering, and dynamical systems theory. In previous works (Celletti and Villa, 1996;Ivanov et al., 1996;Lee et al., 2001;Sarnthein et al., 1998) the same embedding dimension method, has been used for analyzing the electroencephalogram recorded from the brain of different mammals (e.g. rabbits, rats, and humans), and for comparing the results obtained during different experimental (i.e. physiological and pathological, such as epilepsy) conditions. However, the most typical applications of this embedding method are in the area of breathing (Thayer and Moulden, 1997) and heart rhythm analysis (Bettermann and Van Leeuwen, 1998;Yeragani and Radhakrishna, 2003), where the intent of the researches is to identify changes in the electrocardiography recording that could indicate the occurrence of a heart pathology (i.e. arrhythmias). Although there are similar applications in the area of neural function (Jelles et al., 2008;Kannathal et al., 2004;Slutzky et al., 2002;Stam, 2005), we are not aware of studies assessing the dynamical dimension in a specifi c portion of the nervous system.