Abstract
Cullen, Kathleen E. and Daniel Guitton. Analysis of primate IBN spike trains using system identification techniques. I. Relationship to eye movement dynamics during headfixed saccades. J. Neurophysiol. 78: 3259–3282, 1997. The dynamic behavior of primate (Macaca fascicularis) inhibitory burst neurons (IBNs) during headfixed saccades was analyzed by using system identification techniques. Neurons were categorized as IBNs on the basis of their anatomic location as well as by their activity during horizontal headfixed saccadic and smooth pursuit eye movements and vestibular nystagmus. Each IBN's latency or “dynamic lead time” (t_{d} ) was determined by shifting the unit discharge in time until an optimal fit to the firing rate frequency B(t) profile was obtained by using the simple model based on eye movement dynamics,B(t) = r + b _{1} E˙(t); where E˙ is eye velocity. For the population of IBNs, the dynamic estimate of lead time provided a significantly lower value than a method that used the onset of the first spike. We then compared the relative abilities of different eye movement–based models to predict B(t) by using objective optimization algorithms. The most important terms for predicting B(t) were eye velocity gain (b _{1}) and bias terms (r) mentioned above. The contributions of higherorder velocity, acceleration, and/or eye position terms to model fits were found to be negligible. The addition of a pole term [the time derivative of B(t)] in conjunction with an acceleration term significantly improved model fits to IBN spike trains, particularly when the firing rates at the beginning of each saccade [initial conditions (ICs)] were estimated as parameters. Such a model fit the data well (a fit comparable to a linear regression analysis with a R ^{2} value of 0.5, or equivalently, a correlation coefficient of 0.74). A simplified version of this model [B(t) = r _{k} + b _{1} E˙(t)], which did not contain a pole term, but in which the bias term (r _{k}) was estimated separately for each saccade, provided nearly equivalent fits of the data. However, models in which ICs or r _{k}s were estimated separately for each saccade contained too many parameters to be considered as useful models of IBN discharges. We discovered that estimated ICs and r _{k}s were correlated with saccade amplitude for the majority of shortlead IBNs (SLIBNs; 56%) and many longlead IBNs (LLIBNs; 42%). This observation led us to construct a more simple model that included a term that was inversely related to the amplitude of the saccade, in addition to eye velocity and constant bias terms. Such a model better described neuron discharges than more complex models based on a thirdorder nonlinear function of eye velocity. Given the small number of parameters required by this model (only 3) and its ability to fit the data, we suggest that it provides the most valuable description of IBN discharges. This model emphasizes that the IBN discharges are dependent on saccade amplitude and implies further that a mechanism must exist, at the motoneuron (MN) level, to offset the effect of the bias and amplitudedependent terms. In addition, we did not find a significant difference in the variance accounted for by any of the downstream models tested for SLIBNs versus LLIBNs. Therefore we conclude that the eye movement signals encoded dynamically by SLIBNs and LLIBNs are similar in nature. Put another way, SLIBNs are not closer, dynamically, to MNs than LLIBNs.
INTRODUCTION
To generate a saccadic eye movement, motoneurons (MNs) are required to provide signals to the extraocular muscles that overcome the restraining forces of the orbital tissues and drive the eye quickly to its new position and hold it there. Robinson (1964) proposed that the following fourthorder differential equation provides an accurate description of MN discharge in this process
A number of investigators have attempted to simplify Eq. 1 by retaining only the most significant terms. For example, Keller (1973) and Van Gisbergen et al. (1981) proposed a secondorder model including eye position (E), velocity (E˙), and acceleration (Ë) terms. Goldstein and Robinson (1984) and Optican and Miles (1985) also proposed a secondorder form but included as well a pole term [proportional to the derivative of the neuron's firing rate, M˙N(t)] to explain the exponential decay of MN firing following a saccadic eye movement (the “postsaccadic slide”). More recently, Fuchs et al. (1988) demonstrated that a thirdorder model, with only inertia of the globe neglected in the derivation of the model Eq. 1 , provides an adequate description of MN(t). However, Stahl and Simpson (1995) noted that the rejection of the secondorder approximation of Eq. 1 may be premature because of questionable constraints placed on parameters by Fuchs et al. (1988).
The most prominent simplification to Eq. 1
is the firstorder model (Robinson 1970)
For horizontal saccades, excitatory burst neurons (EBNs) project to the ipsilateral abducens nucleus (ABD) to activate it during saccades and act together with a group of inhibitory burst neurons (IBNs) that silence antagonist abducens MNs (Hikosaka and Kawakami 1977; Hikosaka et al. 1978; Igusa et al. 1980; Sasaki and Shimazu 1981; Scudder et al. 1988; Strassman et al. 1986a,b; Yoshida et al. 1982). IBNs and EBNs discharge vigorously or “burst” during ipsilateral saccadic eye movements and are thought to carry similar signals, because IBNs are driven by EBNs (Sasaki and Shimazu 1981; Strassman et al. 1986a) and EBNs appear to be a primary target of the contralateral projections of IBNs (Strassman et al. 1986b).
In recent models of the headfixed saccadic system (Fuchs et al. 1985; Jürgens et al. 1981; Robinson 1975; Scudder 1988), a saccade is controlled by a feedback loop wherein an eye motorerror signal (desired change in eye position minus actual change in eye position) is the “upstream” signal that activates burst neurons (BNs) whereas the “downstream” output of these neurons drives MNs (Fig. 2 B). From the downstream perspective and within the context of the simple model Eq. 2 , BNs provide the eye velocity signal (for review see Fuchs et al. 1985). Accordingly, most prior analyses of BN activity have focused largely on either the relationship between peak eye velocity and peak discharge frequency or global characteristics of the discharge such as the relationship between number of spikes (NOS) and horizontal component of saccade amplitude (Kaneko and Fuchs 1981; Scudder et al. 1988; Strassman et al. 1986a,b; Yoshida et al. 1982). These are static relationships that do not provide information concerning the timevarying relationship between BN spike trains and resultant eye movement trajectories. Only the study by Van Gisbergen et al. (1981) considered timevarying or dynamic relationships. They asked whether the difference in the ondirection (onD) and offdirection (offD) responses of BN signals could provide some of the terms in Eq. 1 other than velocity and examined, with phaseplane trajectories, the influence of acceleration. Using a method that required removing the acceleration term by trial and error, they suggested that the burst discharge is dynamically related to both eye velocity and acceleration.
This paper is the first in a series of three companion papers wherein we analyze in detail the discharges of one class of BNs, IBNs. We have focused on IBNs because 1) their anatomic location is easily defined, which facilitates microelectrode access; 2) they show an excellent static relationship to saccade metrics; and 3) they are driven by EBNs and EBNs appear to be a primary target of the contralateral projections of IBNs (Sasaki and Shimazu 1981; Strassman et al. 1986a,b). In this paper we employ objective computerbased optimization algorithms to characterize the discharges of BNs in the context of Eq. 1 . Our objective is to demonstrate which of the terms in Eq. 1 are important for describing the discharge of this class of premotor neurons. Using algorithms founded in system identification theory (Ljung 1987)we sought to construct the mathematical model that best describes the inputoutput relationship between an IBN's firing rate and properties of the saccadic eye movement trajectory. In the second companion paper of this three paper series (Cullen and Guitton 1997a) we apply these techniques to the analysis of IBN discharges in the headfree monkey to discover whether burst discharges are best correlated to eye, head, or gaze trajectories. In the third companion paper of this series (Cullen and Guitton 1997b) we model IBN firing rate as a function of the upstream motorerror signal (Fig. 2 B) and compare upstream and downstream models in terms of how well they account for the neural discharge pattern.
Our analysis method accomplished the following: 1) provided an objective method for calculating lead time and for comparing the goodnessoffit of different models; 2) provided an indication of whether it is warranted to increase the goodnessoffit by using more complex models that have more free parameters; and 3) permitted a quantitative evaluation of the relationship between estimated initial firing rates and/or biases and movement characteristics such as initial and final eye position, saccade amplitude, and peak eye velocity. The results of our analysis indicate that 1) contrary to the current assumption, primate shortlead IBNs (SLIBNs) and longlead IBNs (LLIBNs) encode eye movement dynamics equally well, and 2) eye velocity–based models that included a bias term (r), which is inversely related to the amplitude of the saccade, provide a good description of BN spike trains.
METHODS
Two monkeys (Macaca fascicularis) were prepared for chronic extracellular recording. The surgical preparation of the animals and the methods used for obtaining the extracellular recordings from brain stem neurons were similar to those previously described by Cullen and McCrea (1993). The IBNs described in this study were identified on the basis of their characteristic burst discharge during rapid orienting horizontal eye movements in the ipsilateral direction and their location in the IBN area caudal and ventral to the ABD (see descriptions for cat in Hikosaka and Kawakami 1977; Kaneko and Fuchs 1981; Yoshida et al. 1982; and monkey in Scudder et al. 1988). Figure 1 illustrates the anatomic distribution of the entire population of IBNs recorded in this study. The location of these neurons was determined by using reconstructions of electrode tracts, stereotaxic locations relative to the ABD, and later histological confirmation. Figure 2 A illustrates the marking lesion made in monkey L that was made just dorsal to the location where two IBNs were recorded. Cells were also characterized as BNs on the basis of their burst responses during quick phases of vestibular nystagmus evoked during horizontal rotation of the vestibular turntable and also by their lack of response during slow phases of nystagmus and smooth pursuit eye movements.
Experimental paradigms and data collection
During the experiments the headfixed monkey was seated in a primate chair that was fixed to the suprastructure of a vestibular turntable. For a juice reward the monkey was trained to orient to a target light that was projected onto a cylindrical screen surrounding the turntable. The spot of light was generated by a HeNe laser that was projected onto the screen via a system of two mirrors, mounted above the monkey's head, whose positions were controlled by galvanometers (General Scanning). The screen onto which the visual target was projected was monochromatic and evenly illuminated. The target was stepped between positions ±5, 10, 20, and 35° horizontal relative to the straight ahead position. The monkeys were also trained to orient to a food target that appeared unexpectedly on any side of an opaque screen facing the monkey (the “barrier paradigm,” see Guitton et al. 1984 for details). This procedure provided us with a rich variety of saccade vectors, starting positions, and, for saccades of the same vector, a variety of velocity profiles.
On average, the saccades in this study were slower than those reported previously in the monkey (e.g., Fuchs 1967). Figure 3 A illustrates the relationship between saccade duration and saccade amplitude for the data set obtained from one animal (monkey H). As expected, saccade duration was related to saccade amplitude; however, there was considerable variation across the data set. For example, the duration of a 30° saccade ranged from 60 to 160 ms. This variability was a result of the wide range of velocities that the monkey utilized to generate a saccade of a given amplitude (Fig. 3 B). The peak eye velocity generated during a 30° saccade ranged from 250 to 800°/s. In general, the saccadic eye velocities in our study were highly variable and often similar to those that were reported in the literature for humans (Robinson 1964). The difference between the metrics of the saccades in our study and other primate studies may have arisen as a result of two factors: 1) the more natural behavioral paradigms employed in our study resulted in monkeys producing less “stereotyped” saccade trajectories to orient to targets (≈70% of the saccades in each neuron's data set) and 2) M. fascicularis were used in this study, whereas most previous studies have used M. mulatta. Although these two types of monkeys have a similar oculomotor range, it is more difficult to obtain consistent behavioral performance in M. fascicularis monkeys (Cullen, unpublished observation).
For this paper BN activity was recorded during rapid voluntary orienting eye movements (saccades) made with the animal's head held fixed. Data for the headfree condition are reported in Cullen and Guiton 1997a. A 386 microcomputer equipped with a data translation DT2825 data acquisition board, a DATAQ waveform scroller board, and a CED1401 peripheral device (Cambridge Electronic Design) was used for controlling visual and vestibular stimuli, for acquiring data, and for online display and analysis of data. Eye movements were recorded by using the magnetic search coil technique. The recorded position signals were lowpass filtered (250 Hz, 8 pole Bessel) and sampled at 1,000 Hz by the computer. The lowpass filtering shifted the recorded position signals backward in time by 2 ms; the resultant delay was taken into account in our subsequent analysis.
Extracellular recordings were obtained with enamel insulated tungsten microelectrodes (7–10 MΩ impedance, Frederick Haer), which were inserted into the cerebellum through a 22gauge guide tube. The guide tube was positioned with an XY stage (Narishige) that was attached to the monkey's recording well. The microelectrode was advanced along the Z axis into the brain stem with a miniature hydraulic drive (Narishige). Spike potentials were amplified and filtered (bandpass 400 Hz to 10 kHz). A unit's isolation was continuously monitored, so that the level of a windowing circuit was set at an appropriate level to generate a pulse coincident with the rising phase of a spike. This pulse was sent to the event channel of the intelligent peripheral (CED1401) device, logging the event times of the neuronal action potentials (see Cullen and McCrea 1993).
Position signals and the single unit activity were simultaneously recorded on Digital Audio Tape (DAT) tape (20 kHz maximum sampling frequency). The DAT recorded data were often replayed offline, so that the unit data could be carefully retriggered [this option was particularly important for the headfree recordings presented in the companion papers (Cullen and Guitton 1997a,b)].
Analysis of BN discharges
During offline analysis sampled eye movement traces were digitally filtered (Matlab) at 125 Hz, because Fourier analysis of ocular saccades has revealed little power above 50 Hz (for review, see Cullen et al. 1994, 1996). Eye velocity was derived digitally from position data. A nonlinear function relating the NOS in the unit discharge to polar angle of the saccade was employed to determine the preferred direction of each IBN (Scudder et al. 1988). The onset and offset of a saccade were defined using a 20°/s eye velocity criterion. Burst duration was defined as the time between the onset and offset of the burst. Burst onset was defined as when the first spike of the burst occurred. Because these neurons frequently generated a few additional spikes following a saccade or gaze shift, burst offset was defined as when 95% of the spikes in the burst had occurred. A spike density function, in which a Gaussian function [5 (SD) ms] was convolved with the spike train, was used to represent the neural discharge. Choosing a Gaussian of this width effectively lowpass filtered the neural discharge so that its frequency content was similar to that of the coincident eye movements (see Cullen et al. 1996 for details). In a previous study (Cullen et al. 1996) we had determined that decreasing the width of the Gaussian by as much as 50% had little effect on the results of our dynamic analysis (latencies and parameters were not markedly changed). Only those saccades in which the amplitude of the vertical position component was less than onethird of the horizontal component were studied; i.e., the amplitude of the horizontal component was between 95 and 100% of the overall movement amplitude. We considered these to be horizontal saccades.
Dynamic models of BN firing rate
We used system identification techniques (Ljung 1987) for the analysis of IBN spike train dynamics. The use of these algorithms permitted optimization across multiple eye trajectories simultaneously and provided a method to objectively determine 1) how well a given model predicts a neuron's discharge pattern; 2) whether increasing the complexity of a model is permissible; 3) what a BN's dynamic latency is; and 4) the relationship between initial conditions (ICs) and quantities describes saccade trajectories. Details of these techniques can be found in Cullen et al. (1996) wherein we describe system identification techniques as they relate to BN analysis. IBN discharges were evaluated in terms of the currently favored schema (Fig. 2 B) (Jürgens et al. 1981) of the saccade control system that states that burst firing rate is related simultaneously to the upstream motorerror signal (Cullen and Guitton 1997b) and to the downstream dynamics of the saccade trajectory (this paper; Cullen and Guitton 1997a). In results of this paper we compare different linear and nonlinear dynamic models that predict IBN firing rate on the basis of downstream signals during headfixed saccades.
METHODS FOR PARAMETER ESTIMATION.
For downstream models, BN discharge was the output and a function of parameters describing saccade trajectory. Therefore these were inverse models because burst frequency was the output of the system, rather than the input as it is in Fig. 2 B. This analysis procedure was necessitated by the estimation techniques that required that the noise be in the output of the model formulation; this is conveniently placed in the BN discharge (Cullen et al. 1996). By using the IBN firing rate as the output, we were able to make comparisons between the fits computed for different plausible models (Tables 2 and 3). For each model optimal fits were made to an ensemble of ∼40 saccades of different amplitudes between 5 and 45° to either the left or right, depending on the neuron's onD. We chose to characterize the discharges of IBNs in terms of horizontal eye movements, because IBNs project to the MNs that innervate the muscles that are responsible for producing horizontal eye movements (i.e., the lateral and medial recti). Note, in a preliminary analysis we analyzed the responses of several cells along their optimal response direction (see Results) and found that goodnessoffit, relative importance of estimated parameters, and the optimal latencies estimated in our models were not markedly changed. Generally, the starting positions of the saccades in each data set were well distributed over a range of ±20° on either side of the center position. Optimal fits to the neuron's discharge during offD saccades were carried out separately.
ESTIMATION OF DYNAMIC LATENCY.
We shifted the burst in time relative to the saccade by the latency (the optimal dynamic delay t_{d} ) that provided the best fit of the BN firing rate when a simple but accurate downstream model was used (see section on lead time in results). It is important to note that the entire portion of the temporally shifted spike train that was coextensive with the saccade duration was used to fit the model. Thus every spike within this interval contributed to the latency calculation. When the BN was shifted by the appropriate t_{d} , the main portion of the burst was typically aligned with the duration of the saccade. Consequently, the action potentials discharged by the neuron before the main portion of the burst (the burst “prelude”) were typically rejected by the optimization algorithms and were not included in the fit. The estimated delay (t_{d} ) was then applied to the other more complex models (i.e., it was a fixed parameter) (see Cullen et al. 1996 for details).
COMPARISON OF MODEL FITS.
The goodnessoffit of each model was determined by considering the root mean square (RMS) of the fit errors and the variance accounted for (VAF). The VAF was defined as [1 − (var/SD)] where var is the variance of the fit's error and SD is the standard deviation of the data set about its mean. Using VAF allowed us to normalize the goodness of a model's fit across neurons. It should be noted that the VAF for a linear model is equivalent to the square of the correlation coefficient (R; i.e., a VAF of 0.5 indicates that 50% of the variability in the unit's firing rate is explained by the model and would correspond to R = 0.71). We tried models of different complexities. A more complex model has more adjustable parameters and hence has the potential for providing a better fit to the data. The question then arises as to whether increasing model complexity is warranted. As described in Cullen et al. (1996) we answered this question objectively by using the Bayesian information criterion (BIC) (Schwarz 1978); a “cost index” that indicates whether increasing the complexity of a model is justified when the accompanying decrease in the error of the fit is taken into account. We computed the BIC for each model estimation; a decrease in the BIC value indicated that an increase in model complexity was warranted. If this index did not decrease when a more complex model was used, then the model was considered to be no better at describing the data than the more simple model.
RESULTS
We present data from twentyeight IBNs that were recorded during saccadic eye movements in two headfixed monkeys. A characterization of the activity of the cells in terms of the trajectory of headfixed saccades is presented in this paper. In our companion papers we 1) analyze the activity of these neurons in the headfree condition and compare the results with the headfixed data (Cullen and Guitton 1997a) and 2) consider IBN discharges in terms of their response to a motorerror signal (Cullen and Guitton 1997b).
General discharge characteristics
Neurons were identified as IBNs on the basis of their anatomic location (see Methods) and activity during saccadic and horizontal smooth pursuit eye movements and horizontal vestibular nystagmus. The neurons were identified on the basis of their characteristic burst discharge during horizontal saccadic eye movements in the ipsilateral direction. We determined that the cells were not the “burstdriver” neurons that were reported in a similar region in the cat (Ohki et al. 1988) by verifying that they were inactive during slow phases of vestibular nystagmus in both directions. It was also verified that these neurons were inactive during smooth pursuit eye movements and therefore that none could be categorized as the reticular formation smooth pursuit cells previously described by Tomlinson and Brance (1991). The neurons in this study were categorized as either SLIBNs (n = 16) or LLIBNs (n = 12), depending on whether the mean period between the onset of the first spike and the onset of eye velocity was ≤15 ms or >15 ms, respectively; we used this arbitrary definition to be compatible with the prior IBN study of Scudder et al. (1988).
The discharge of a typical SLIBN (H0415) is illustrated in Fig. 4, A–C. This neuron generated a strong burst of firing in the onD (ipsilateral) and weak burst in the offD (contralateral), with the first spike leading ipsilateral saccadic eye movements on average by 11 ± 2.8 (SD) ms (Fig. 4 A). In Figs. 4 and 5 and all subsequent figures showing spike frequency, the spike discharge was shifted to the right by a time equal to the calculated dynamic lead time of each neuron (see Estimation of lead times). The duration of this neuron's burst was tightly correlated with the duration of the saccade (slope = 0.9; R = 0.91) and the total NOS during a burst was proportional to the horizontal component of the amplitude of the saccadic eye movement (slope = 1.5; R = 0.92). During smooth pursuit eye movements this cell only discharged during saccades and did not discharge during slowphase smooth pursuit movements in either the ipsilateral or contralateral directions (Fig. 4 B). During vestibular stimulation this cell discharged only during quick phases and did not discharge during ipsilateral or contralateral slow phases (Fig. 4 C).
Figure 4, D–F, illustrates the discharge of a typical LLIBN (H0421). This neuron generated a preamble of lower frequency discharge followed by a burst of firing, with the first spike in the preamble leading ipsilateral saccadic eye movements by 21 ± 7.9 ms (Fig. 4 D). The duration ofthis neuron's spike train was correlated with the duration of the saccade (slope = 0.9; R = 0.74) and the total NOS generated during the train was proportional to the horizontal component of the amplitude of the saccadic eye movement (slope = 0.9; R = 0.62). As was the case for SLIBNs, this neuron did not discharge in response to slowphase smooth pursuit or vestibularly induced eye movements in either the ipsilateral or contralateral directions (Fig. 4, E and F).
IBN preferred direction
The IBNs in this study burst vigorously for ipsilaterally directed saccades, and they typically discharged far fewer spikes during contralateral (offD) and pure vertical saccades. Figure 5 shows examples of discharges of three neurons we selected to illustrate our data throughout this paper and the companion papers (Cullen and Guitton 1997a,b): SLIBN L0702, LLIBN H0409, and LLIBN H0925. The latter neuron had the best offD discharge in our population. Note the clear correlation between burst frequency and saccade velocity for discharges in the onDs and offDs for H0925 and in the onD for cells L0702 and H0409. A precise measure of a neuron's optimal saccade direction was calculated by fitting a nonlinear function to relate the NOS in a unit's discharge to the polar angle of the saccade (see Methods). The saccade amplitudes used to determine the preferred direction of each cell in the study ranged from 15 to 25°. The preferred directions in our population (Fig. 6) varied from +29° (upward) to −30° (downward). The mean preferred direction was −1.0 ± 15.2° closely aligned with the horizontal stereotaxic plane (∼15° upward from the plane of the horizontal canals). Our analysis of eye movements that were horizontal, or nearly so, is consistent with the mean preferred horizontal direction of IBNs (see Methods).
Estimation of lead times
The time by which a BN's discharge led eye saccades was calculated by two methods. Histograms illustrating the results of applying each method are shown in Fig. 7. In Fig. 7
A the lead time was determined by calculating the period between the onset of the first spike and the onset of eye velocity. This method was used in many previous studies of IBNs and as a criterion to categorize the cells as either SLIBNs or LLIBNs (e.g., Scudder et al. 1988). In the second method (Fig. 7
B), by using our optimization algorithms described in detail in Cullen et al. (1996), a neuron's dynamic lead time (t_{d}
) was left as an unknown and determined by shifting progressively the unit discharge in time relative to the saccade until a best fit was obtained for models 2d and 6d (Table 2)
The relationship between the lead times calculated by using each method is plotted in Fig. 8. The values obtained from the two methods were not well correlated for the 28 primate neurons we analyzed (R = 0.51). However this relationship improved once a single outlier (the value within parentheses) was removed from the analysis (R = 0.67). A continuity can be seen between the dynamic lead times estimated for SLIBNs (♦) and LLIBNs (⋄). The timeoffirstspike method yielded a mean lead time for SLIBNs that was statistically similar to the one using dynamic analysis with Eq. 3 (11.6 ± 3.8 ms vs. 11.8 ± 2.7 ms, respectively). However, for LLIBNs the firstspike method yielded significantly longer lead times (24.2 ± 9.7 ms vs. 15.5 ± 2.2 ms, respectively). For the entire population of neurons, including SLIBNs and LLIBNs, the mean lead times were 17.0 ± 9.3 ms and 13.4 ± 3.1 ms for the firstspike and dynamic methods, respectively. The difference between these values was significantly different (Student's ttest,P < 0.01). It is noteworthy that the mean lead time estimated dynamically for neurons classified as SLIBNs was significantly less (P < 0.001) than that estimated, using the same method, for neurons classified as LLIBNs [11.8 ± 2.7 ms vs. 15.5 ± 2.2 ms, respectively, a difference of 3.7 ms to be compared with the 8 ms estimated by Scudder et al. (1988)]. Only five neurons (4 LLIBNs and 1 SLIBN) had dynamic lead times longer than 15 ms.
Relationships between IBN activity and saccade metrics
Linear relationships between the duration of IBN discharges and ipsilaterally directed saccadic eye movements were reported in a number of previous studies and our data are comparable. Figure 9, A and B, illustrates the relationship between burst duration and saccade duration for neurons classified as SLIBNs and LLIBNs, respectively. The insets show data obtained from neurons that we selected as examples and whose data are illustrated throughout the paper. The heavy lines show the population means whose slopes, intercepts, and correlation coefficients are given in Table 1. The slopes for SLIBNs and LLIBNs were not significantly different (P = 0.11), but the intercepts were (P < 0.05). In Table 1 we also show data from Scudder et al. (1988). Both data sets concur; the average burst duration of SLIBNs closely approximates the duration of the associated saccade, whereas for LLIBNs the discharge preamble produces spike trains that are longer than saccades. The present analysis suggests that the duration of this preamble, for a given cell, is constant (because the slope is ∼1) and independent of saccade duration (and amplitude).
A number of studies have reported a strong relationship between the NOS in an IBN's discharge and the horizontal amplitude of ipsilaterally directed saccades. To obtain the NOS in the portion of the burst that directly drives saccades, we shifted the spike train in time by the dynamic lead time (t_{d} , calculated by using model 2d of Table 2) and counted the NOS during the discharge period that was aligned with each saccade (e.g., between the vertical dashed lines in Fig. 4 A). The NOS was proportional to the amplitude of the horizontal component of the saccadic eye movement for each neuron. (Because the saccades we analyzed were nearly horizontal, using the overall amplitude did not influence this conclusion.) This relationship is illustrated for SLIBNs and LLIBNs in Fig. 10, A and B, respectively. The heavy lines show the population means whose slopes, intercepts, and correlation coefficients are shown in Table 1, as well as the equations given by Scudder et al. (1988). The slope and intercept of the mean lines were not significantly different for SLIBNs and LLIBNs (P = 0.41 and 0.11, respectively) and are also similar to Scudder et al. (1988). We conclude that our method of defining the interval over which NOS were counted, by using an objective estimate of lead time, preserves the strong relationship that was observed in previous studies between total NOS in a burst and the horizontal amplitude of a saccade.
Previous studies also reported a relationship between IBN peak firing frequency and peak saccade velocity (Scudder et al. 1988; Yoshida et al. 1982). We did not find this relationship for SLIBNs and LLIBNs (Table 1; mean R = 0.39) to be as robust as that in Scudder et al. (1988); it was significant for 38 and 8% of our SLIBNs and LLIBNs, respectively. This difference may be due to the lessstereotyped saccade profiles that were elicited during some of our protocols (see Methods and discussion).
Dynamic models of IBN discharges during headfixed saccades in OND
The foregoing analysis of IBN lead times and relationships with saccade metrics has 1) characterized our cells as typical IBNs by using criteria similar to those of other studies and 2) introduced and employed a dynamic estimate of discharge frequency lead time, which is further utilized in the following text to evaluate different models that relate IBN spike train dynamics to eye movement parameters during onD saccades. Note that the action potentials discharged by the neuron before the main portion of the burst (the burst prelude) and after (the tail) were typically excluded by the optimization algorithms (see Methods). The resulting fits are illustrated throughout the paper for the example SLIBN (L0702) and LLIBN (H0409).
Table 2 shows the eight downstream models we investigated. It also provides the means and SDs of the parameters estimated for each model for SLIBNs, LLIBNs, and the entire population of IBNs. Table 3 gives the VAF and BIC produced by each of these models for each of the 28 neurons in this study as well as the means and SDs for SLIBNs, LLIBNs, and for the entire population.
We will now consider each model in turn beginning with the simplest.
BURST FREQUENCY PROPORTIONAL TO EYE VELOCITY.
Model 1d in Table 2 is the simplest. This often used model addresses the question of whether instantaneous burst frequency is proportional to eye velocity. It is based on the wellknown relationship (discussed previously) between peak discharge frequency and peak eye velocity (Kaneko and Fuchs 1981; Scudder et al. 1988; Strassman et al. 1986a,b; Yoshida et al. 1982).
The VAF and BIC values for model 1d are shown in Table 3. This model was a poor predictor of the data as indicated by the low VAFs. Indeed, a negative VAF in 12 of 28 of the neurons implies that the variance of this model's fit error was actually larger than the variance of the data about a constant mean level.
ADDITION OF BIAS TERM.
Model 2d is an extended form of 1d and can be viewed as a firstorder linear approximation to a nonlinear relationship between B(t) and E˙ where r is the bias term.
When this bias term was added to the model (Table 2), the VAF for the entire population increased very significantly from 0.04 to 0.30 (Table 3). Note that this model actually provided a fit to the data that was on average as good as that described by a correlation coefficient of 0.55 in a linear regression analysis. The accompanying decrease in the average BIC value (from 10.23 to 9.60) indicates that the addition of a bias term was warranted.
Averaged measures of goodnessoffit were also calculated separately for SLIBNs and LLIBNs for the models 1d [Table 3: VAF (shortlead) = 0.05; VAF (longlead) = 0.03] and 2d [Table 3: VAF (shortlead) = 0.32; VAF (longlead) = 0.28]. There was no significant difference in the VAF by each of these models nor by any other model we tested for the two subpopulations of cells; accordingly in Table 2 we pooled the results for both SLIBNs and LLIBNS and display the average value for each parameter in each model. The considerable improvement in the fit of IBN firing rate that occurred with the addition of the bias term is illustrated for three example saccades for SLIBN L0702 (Fig. 11 A) and LLIBN H0409 (Fig. 11 B). The goodnessoffit to the SLIBN is particularly striking.
IMPORTANCE OF ACCELERATION TERM.
The next downstream model we investigated utilized a simple approximation to the hypothesis of Van Gisbergen et al. (1981), that unit firing can be represented by a nonlinear function of eye velocity plus an eye acceleration term. As a first approximation to this proposal, the nonlinearity in velocity was (as in Eq. 3 ) simply represented by the bias term (r; Table 2: model 3d). For all our neurons the addition of the acceleration term (b _{2}) had little effect on increasing the VAF (or decreasing the BIC); the addition of an acceleration term increased the mean VAF from 0.30 to 0.31 (Table 3). The same was true when the SLIBNs and LLIBNs were considered separately. Examples of the marginal improvement in the fit of the IBN firing rate that accompanied the addition of an acceleration term are illustrated in Fig. 11.
NONLINEAR VELOCITY TERMS.
A more complex nonlinear function of eye velocity was tested that included a thirdorder nonlinearity in addition to an eye acceleration (Table 2: model 4d). The addition of the higherorder eye velocity terms (d _{1} and d _{2}) decreased the error of the model slightly (Table 3: mean population VAF = 0.33). Taken together the results so far show that in downstream models a bias term is much more significant than either acceleration or nonlinear velocity terms for representing IBN discharges.
ESTIMATION OF POLE TERM.
A model was also investigated which in addition to eye velocity, acceleration and bias terms allowed for the estimation of a pole term in the system (the derivative of the IBN firing rate; Table 2: models 5d and 6d). This corresponds to the slide component of the MN's response described in Eq. 1 by the term cM˙N(t). The addition of this term provided a slightly better model fit than model 3d when the ICs in firing rate (necessitated by the pole term) were taken from the data (Table 3: VAF = 0.31 in model 3d vs. 0.37 in model 5d). However, when the ICs were estimated as parameters (Table 2: model 6d), the fit improved considerably (Table 3: VAF = 0.37 in 5d vs. 0.54 in 6d). For model 6d the fits to the data were as good, on average, as those described in a linear regression analysis by R values of 0.74 and 0.73 for SLIBNs (VAF = 0.55) and LLIBNs (VAF = 0.53), respectively. Examples of fits to firing rate are shown in Fig. 12 for the same example neurons and saccades as in Fig. 11. The improvement in fit when ICs were estimated as parameters (2nd panel from top) rather than taken from the data set (top panel) can be appreciated by visual inspection. In Fig. 12 the numbers above each firing frequency profile give the initial frequency B(t _{0}).
A model of the form 6d requires the estimation of many more parameters (equivalent to the number of saccades included in the analysis) than one in which the ICs were taken directly from the data (model 5d). To determine whether increasing the model complexity to this extent (from 4 to 44 parameters) was warranted, it was necessary to consider the BIC. Comparison of the BIC values for models 3d, 5d, and 6d (Table 3) shows that the addition of ICs as parameters to be evaluated was warranted (mean BIC = 9.57, 9.41, and 8.84, respectively). Model 6d, however, suffered from unrealistic values of the bias term and this will be considered in the discussion.
ESTIMATION OF VARIABLE BIAS.
In model 6d, the average time constant of the decaying ICs is long relative to saccade duration and given by the value c = 2.4 s (see Table 2; discussion). Thus the ICs were effectively constants where the model estimated a different constant value for each saccade that was the sum of the estimated value of the bias and the constant IC. This observation led us to specifically test a model in which only one bias term was estimated separately for each saccade (Table 2: model 7d). This variable bias model had nearly as many parameters as did model 6d (42 vs. 44, respectively); however, the method of estimation was simpler for a model of this form than for models that include a pole term (see Cullen et al. 1996). Example fits obtained using model 7d are illustrated in Fig. 12 (3rd panel from top). This model provided fits that were nearly as good as those given by model 6d (Fig. 12; Table 3: VAF = 0.54 for model 6d vs. 0.52 for model 7d). The possible physiological implications of the variable bias are considered in the following section and in discussion.
ICs AND BIASES DEPEND ON SACCADE PARAMETERS.
For 14 cells—the majority of SLIBNs (9/16) and many LLIBNs (5/12)—the values estimated for ICs in model 6d were inversely and significantly correlated with the metrics of saccades. This was not the case for the example neurons shown in Figs. 11 and 12, and so these relationships are illustrated for another LLIBN (H0925) in Fig. 13, A and B. The estimated ICs for this cell were well correlated with saccade amplitude (R = −0.86) and less strongly related to peak eye velocity (R = −0.62). For those neurons with significant correlations, ICs were better correlated with saccade amplitude (mean SLIBNs, R = −0.70; LLIBNs, R = −0.67) than with peak saccade velocity (mean SLIBNs,R = −0.63; LLIBNs, R = −0.59). Note that in all these neurons the value of the ICs varied inversely with amplitude and peak velocity (as in Fig. 13).
For 9 of the 14 cells whose ICs were correlated to metrics, the bias term estimated for model 7d (Fig. 13, C and D) was also linearly correlated with the amplitude and/or peak velocity of the saccade. The estimated biases for cell H0925 were well correlated with saccade amplitude (R = −0.89) and less strongly related to peak eye velocity (R = −0.56). The bias term was largest for small saccades. This implies that model 2d, which has a fixed bias, averaged across all saccade amplitudes should be best at predicting mediumamplitude saccades. Furthermore, model 2d should underevaluate the discharge frequency of small saccades and overevaluate that of large ones. This is precisely what Fig. 14 shows for neuron H0925.
We also investigated whether a relationship existed among either the estimated biases and/or ICs and the position of the eye before and/or following a saccade. For two SLIBNs there was a significant relationship among the estimated ICs, as well as the estimated biases, and the position of the eye before the onset of the saccade; however, no relationship was observed with final position for any IBN in this study. No relationships were ever observed among the ICs when they were taken directly from the data (model 5d) and any of the measured saccade parameters (amplitude, peak velocity, or initial or final eye position).
EFFECT OF SACCADE AMPLITUDE ON IBN SPIKE TRAIN DYNAMICS.
Our finding that estimated ICs and biases were often correlated with saccade amplitude led us to construct a more simple model that included an amplitudedependent term. This model (Table 2: model 8d) is a simple extension of model 2d, to which a term that is proportional to saccade amplitude (r _{1}ΔE) was added. The addition of the eye amplitude term to model 2d resulted in comparable improvements in VAFs of LLIBNs (Table 3: VAF = 0.28 in model 2d vs. 0.32 in model 8d) and SLIBNs (Table 3: VAF = 0.32 in 2d vs. 0.36 in 8d). On average, model 8d was slightly better at describing IBN discharges than the nonlinear function of eye velocity in model 4d (Table 3: VAF = 0.34 in 8d vs. 0.33 in 4d). However, model 8d did not fare as well as 7d, but it had 38 less parameters.
EFFECT OF EYE POSITION.
The final two models that were tested addressed the issue of whether eye position was relevant in estimating IBN discharges. The simplest model included two terms, an eye position term and a bias. For each neuron this model was far less effective at describing the discharges of cells than the simple model of Eq. 3 that included an eye velocity term and a bias (VAF = 0.05 vs. 0.30 in model 2d). Furthermore, the addition of a position term to model 2d resulted in only a negligible increase in the VAF (VAF = 0.31).
Downstream model errors for SLIBNs vs. LLIBNs
We have shown that, as a population, SLIBNs have a shorter dynamic lead time than LLIBNs. This is compatible with a linear signal processing stream wherein it is possible that LLIBNs project to SLIBNs that in turn project to MNs (Fuchs et al. 1985). However, we argue here and in the companion paper (Cullen and Guitton 1997b) that this is not the case. To determine whether SLIBNs encode downstreameye movement dynamics better than LLIBNs, we examined the relationship between lead time and the VAF values estimated for each downstream model. The VAF was not lower for those neurons that had shorter lead times; hence, SLIBNs do not appear closer than LLIBNs to the motor output for any of the eight models that we tested. The relationship between lead time and VAF for model 8d is illustrated in Fig. 15.
OFFD discharges
The offD discharge of IBNs is of considerable interest because it will subtract from the excitatory drive provided by EBNs on agonist MNs. It is therefore important to provide a quantitative evaluation of this effect if we are to understand how MN signals are generated (Van Gisbergen et al. 1981).
Figure 5 shows the offD of the three neurons that we used extensively as examples in this paper. The offD of SLIBN L0702 and LLIBN H0409 consisted of an occasional one or two spikes and could not be analyzed. It is important to note that 50% of each of our SLIBN and LLIBN populations gave similar weak offD discharges. The other 50% provided offD discharges that varied in the quality of their relationships to eye dynamics; of these, cell LLIBN H0925 exhibited the best model fits, as can be appreciated from Fig. 5. We will consider this neuron in detail.
For cell H0925, the dynamic lead time in the offD was 11 ms compared with 17 ms in the onD. The onD and offD data sets contained a similar distribution of saccade sizes. Consequently, the difference in onD and offD dynamic lead times does not appear to be the result of a difference in the amplitude of the saccades that were included in the analysis. As for the correlation coefficient in the linear relation between NOS and amplitude, the values were 0.86 and −0.71 in the onDs and offDs, respectively. Table 4 compares for this neuron the characteristics of the model fits for each model we have tested. The VAFs in the onD and offD were generally comparable. The comments made in an earlier section on the relative accuracy of each model, as model characteristics and complexity were changed, apply to the offD. Thus models 6d and 7d best describe the offD response of cell H0925. Indeed, in some cases the offD fits were superior to the onD fits, as exemplified by model 6d where the VAF was 0.57 for the offD compared with 0.50 for the onD. The offD fit in this case is equivalent to a linear regression coefficient (R) of ∼0.76.
To compare parameters in the onD and offD models we will focus on models 2d, 7d, and 8d because of their similarity in structure and also because we later propose that 8d is the simplest, most accurate and physiologically plausible representation of the cell's onD discharges (see Discussion). In model 2d the absolute value of the velocity gain, b _{1} = −0.83, was ∼10% lower in the offD compared with the onD and, most strikingly, the bias was far less and essentially zero in the offD. This latter point is emphasized by the same VAF for the offD in models 1d and 2d. Similarly in model 8d, the absolute value of the velocity gain was ∼20% lower in the offD (b _{1} = −0.86) compared with the onD, and the bias was nearly zero in the offD. The eye amplitude dependent gain factor (r _{1}) had considerably less effect on the model fit in the offD (VAF = 0.33 in model 8d vs. 0.32 in model 2d) than in the onD (VAF = 0.45 in 8d vs. 0.32 in 2d). Accordingly, the amplitude gain r _{1} was considerably less in the offD (1.5) than in the onD (−7.5). In model 7d the velocity gain b _{1} = 0.90 was also ∼10% lower in the offD, and unlike in the onD for this cell, no relation was found between bias values and saccade amplitude. This result is compatible with the relatively small value of r _{1} estimated for model 8d in the offD. Although the average bias in 7d was equal to that in 2d (i.e., essentially 0), the fact that the bias was allowed to vary slightly improved the VAF as noted above.
DISCUSSION
Reproducibility of analysis
The main objective of this paper was to provide an accurate description of IBN burst discharges to better understand the signal carried by these cells and how it relates to MN discharge during saccade generation. Our analysis, by using system identification techniques, has provided an objective method for comparing the goodnessoffit of different downstream models, as well as an indication as to whether increasing the complexity of a model by adding parameters is warranted.
Before going on to discuss the significance and implication of our results, it is important to demonstrate that a given model fit can be extended beyond the data set for which it was constructed. This is verified in Fig. 16 for models 8d and 7d applied to SLIBN H0923. Figure 16 A shows three examples of each of the best and worst fits applied to the original data set. The average VAF by the three best fits was 0.50 and 0.60 in models 8d and 7d, respectively, compared with 0.28 and 0.42, respectively, for the three worst fits. Figure 16 B shows the fits to a new data set from the same models as used in Fig. 16 A and with the same parameters (r _{0}, r _{1}, and b _{1} in model 8d and b _{1} in model 7d). In the latter model the biases were fit as parameters to the new set of 40 saccades and are independent of those derived in the original data set. The VAF by models 8d and 7d for this new population was 0.35 and 0.50, respectively, compared with 0.37 and 0.59 in the original data set. The average VAF by the three best fits was 0.51 and 0.59 in models 8d and 7d compared with 0.11 and 0.32, respectively, for the worst fits. These data show that the models predict IBN firing in the new data set almost as well as in the set used to calculate the parameters.
Comparison of results of our metric analysis with other studies
Our analysis of the relationship between saccade metrics and IBN discharges revealed that the cells in our sample were generally similar to the primate IBNs described in the study of Scudder et al. (1988). Each neuron demonstrated a monotonic increase in the NOS in a burst (excluding preludes) with increasing amplitude of the horizontal component of a saccade [R = 0.79 this study vs. 0.85 Scudder et al. (1988)]. Consequently, in the terminology of Hepp and Henn (1983), our SLIBNs as well as LLIBNs would have been classified as “directional” bursters rather than “vector” bursters. In addition, we found that the preferred direction of our population of IBNs was closely aligned with horizontal (onD was −1.0 ± 30°). This result is in agreement with that of Scudder and colleagues (1988) who reported a similar spread in the ranges of preferred directions for their sample of LLIBNs and SLIBNs. In both studies, the range of preferred directions is much smaller than that previously reported in squirrel monkeys (Strassman et al. 1986b).
Several previous studies have reported a strong relationship between IBN peak firing rate and peak eye velocity during saccades. Our use of a spike density function rather than reciprocal intervals to represent unit firing rate resulted in a filtering of the firing frequency profile; we expected to find even better correlations. We were surprised to find a poor relationship between IBN peak firing frequency and saccade peak velocity (mean R = 0.39). This difference could stem from the relatively nonstereotyped skewed saccade profiles (e.g., left saccade in bottom panel of Fig. 11 B; right saccade in bottom right panel of Fig. 16 B) elicited during one of our protocols (the barrier paradigm; see Methods) where for a given amplitude the peak velocity was quite variable.
Significance of dynamic lead time
Previous studies of BNs determined lead times based on a variety of methods such as the time of first spike or transition from low to high frequency (see a review in Hepp et al. 1989). One advantage of the present study was that an objective dynamic lead time was calculated for each unit, which accounted for the average influence of each spike on the saccade trajectory. This dynamic estimate of lead time was then implemented in our analysis of different models to correct for latency. A neuron's dynamic lead time was determined by shifting the unit discharge in time (t_{d} ) until an optimal model fit was obtained for the simple but accurate dynamic model of Eq. 3 : B(t) = r + b _{1} E˙(t − t _{d}). For LLIBNs the dynamic estimate provided a significantly lower value for lead time than that obtained by measuring the onset of the first spike. For SLIBNs, lead times determined by the two methods were the same. As seen in Fig. 8, the lead times calculated by using each method were moderately well correlated. Interestingly, neurons arbitrarily classified as SLIBNs on the basis of the timing of the onset of the first spike still had significantly smaller (by 3.7 ms) dynamic lead times than neurons classified by using the same criterion as LLIBNs. We will consider this in a subsequent section. As expected given this small 3.7ms difference, there was a considerable amount of overlap among the dynamic lead times estimated for SLIBNs and LLIBNs (see Fig. 8); 8 of 12 neurons classified as LLIBNs with the firstspike method were now in the shortlead range.
Because of the mechanisms leading to cell depolarization the initial few spikes of an excitatory burst impinging on MNs should have a different influence than those spikes occurring in the middle of the burst. This phenomenon is also dependent on initial firing rate. In past studies the objective duration of the LLIBN prelude has been difficult to evaluate. Our analyses reveal that on average the time difference between the firstspike and the dynamic lead time was ∼9 ms. Put another way, for LLIBNs 9 ms of activity does not carry information related to the dynamics of saccades, but is rather a preparatory or anticipatory response. As concerns the influence of spikes well within the burst, it is theoretically possible that t_{d} can vary throughout the saccade. We have not tested for this.
Evaluation of models and their estimated parameters
Table 2 gives mean values and SDs of the parameters that were estimated in each of the nine downstream models tested for SLIBNs and LLIBNs and for the entire population taken together. Specific values for our example neurons are provided in the model fits illustrated in Figs. 11 and 12. We have shown that the most influential terms for predicting IBN firing rate are eye velocity and a bias. A model (2d) comprised of these two terms provided a fit that corresponded to a correlation coefficient (R) of 0.55 in a linear regression analysis (VAF = 0.30). The sign of the estimated eye velocity term was positive (ipsilateral) for every neuron in our analysis, compatible with previous studies demonstrating increases in IBN firing rate accompanying increases in eye velocity in a direction ipsilateral to the neuron's location. The estimated bias term was positive in all but the model of Eq. 4 (model 6d) where ICs were estimated as parameters; this result will be considered later.
For the population of neurons, the estimated bias and eye velocity gain terms did not change significantly between models 2d and 3d. This is undoubtedly due to the small acceleration term in model 3d. Indeed, the mean estimated value of b _{2} in model 3d was very low (−0.0004 ± 0.0009) and was unexpectedly in the off, or contralateral, direction. We have shown earlier (Cullen et al. 1996) that the magnitude of the acceleration term is critically dependent on the choice of lead time. For example a smaller than optimal lead time predicts a larger acceleration term because more of the initial burst frequency increase must be accounted for by early changes in eye trajectory. Our objective estimate of best lead time is responsible for minimizing b _{2}. The addition of a position term to model 2d only very slightly improved the fit of IBN discharges. The mean estimated eye position gain was low (−0.04 ± 1.3) and also was in the contralateral direction.
The existence of a bias term around which IBN discharges are modulated during saccades was included in previous models of the brain stem burst generator (Scudder 1988; Van Gisbergen et al. 1981). The average value of the bias term in our estimations using model 2d was 223 spikes/s; relatively close to the value in the Scudder 1988 model (160 spikes/s).
The bias value was significantly lower in the fourthorder nonlinear representation of the Van Gisbergen et al. (1981) model (model 4d) than in the linearized version (model 3d) and the eye velocity gain term b _{1} was correspondingly higher in 4d than in 3d. These changes are compatible with an improved fit by the nonlinear model; but the improved VAF by model 4d with respect to 3d of ∼7% (Table 3) is less than we expected based on the nonlinear relationships between B(t) and E˙(t) shown by Van Gisbergen et al. (1981). In fact model 8d, a simple extension of model 2d, which included an amplitudedependent term, better described our IBN discharges than the nonlinear function of eye velocity in 4d. A number of these points are illustrated in Fig. 17, which shows phaseplane plots [B(t) vs. E˙(t − t _{d})] for cell H0925. In Fig. 17 A the dynamic latency was chosen as the optimal one (17 ms) for this cell and the phaseplane plot trajectories of each saccade are close, on average, to straight lines. Figure 17 B shows the phaseplane plots, based on optimal latency, of three saccades of different amplitudes, all of which were included in Fig. 17 A. Model 2d (——) approximates well the average phaseplane characteristics, but one can see from Fig. 17 B that variations among saccades that are dependent on saccade amplitude could be fit better by including a saccade amplitudedependent term in the dynamic model (i.e., model 8d) or by allowing for the estimation of variable biases (i.e., model 7d). The nonlinear model (model 4d) for this cell also is plotted in Fig. 17 A (– – –) and illustrates well the small nonlinear contributions that we have observed for all cells over the velocity range studied.
When a pole term was added to model 3d to achieve model 5d, the bias increased significantly (ICs were taken from the data). However, the improved fit with 5d was small compared with that offered by 6d when ICs were fit as parameters. Model 6d was the most complex model tested and included velocity, acceleration, and a pole term; the ICs were estimated as parameters. For this model the bias value varied greatly between neurons and the mean was an extremely low value (Table 2: −1,795 ± 2,058). The values of the other parameters in this model were also extremely variable. A large negative bias resulted from the algorithm seeking to offset the large IC transient that is itself a bias, given that the average time constant of the ICs is long relative to saccade duration and given by the value c = 2.4 s. Thus negative biases, in this context, do not have an actual physiological significance; they add to the relatively constant (during a saccade) positive IC transient to produce an overall bias close to that given by model 7d.
To pursue these arguments, note in model 6d that c ≈ b _{2}. The significance of this can be appreciated by noting from model 2d that B˙(t) ∝ Ë(t − t _{d}), approximately, and substituting cB˙(t) = b _{2} Ë(t − t _{d}) in 6d yields model 7d. Taken together, these arguments indicate that in 6d the bestfit algorithm has selected 1) c to null the acceleration term and 2) the bias to null the effect of c. Thus model 6d essentially reduces to 7d (which has neither acceleration nor pole terms) and that provides an excellent description of IBN firing rate. This further substantiates the validity of our algorithms and analysis procedure because 7d, with no pole term, and 6d, with a pole term, were analyzed in quite different ways (Cullen et al. 1996). Model 7d has no pole term, is of comparable complexity to 6d, and is easier to estimate. This formulation provided more meaningful estimates of bias (always positive). However, models 6d and 7d contained too many parameters to be considered as useful models of IBN discharges. Therefore the correlation, which we frequently observed between estimated ICs and biases, led us to construct a more simple model that included an amplitudedependent term (model 8d). This model actually described IBN spike trains slightly better than a thirdorder nonlinear function of eye velocity (model 4d). The sign of the amplitudedependent term in this model was negative for 26 of 28 neurons; that is, this component of the firing rate was inversely correlated with saccade amplitude. Given the small number of parameters in model 8d (only 3) and its ability to fit the data, we suggest that it provides the most valuable description of IBN discharges during saccades.
The amplitudedependent component of the IBN burst that was revealed by the r _{1} term in model 8d and by the amplitudedependent biases estimated in model 7d was not expected. As illustrated in Fig. 17 B, this relationship does not result from piecewise linear fits to the nonlinear relation between B and E˙ (where the bias is the intercept) because the nonlinear terms are very small. Both Van Gisbergen et al. (1981) and Scudder et al. (1988) found that more spikes are needed for a given amplitude change in small (<2°) saccades than for an equivalent amplitude change in large saccades. This result is compatible with a nonlinear relation between B(t) and E˙(t), which states that the effective viscosity is higher for small saccades. However, these arguments do not apply to our data set because we have analyzed saccades larger than 5°; Fig. 17 B shows that there are still differences for cell H0925 between phase plane plots for 14, 21, and 36° saccades. We argue in companion paper III (Cullen and Guitton 1997b) that the amplitudedependent component of the IBN burst may be due to the influence of a signal from the superior colliculus on IBNs.
Acceleration term does not improve fits of IBN discharges
Our finding that the addition of an eye acceleration term did not improve the fit to IBN discharges was unexpected in light of the previous results of Van Gisbergen et al. (1981). In their analysis, these investigators correlated movement trajectory with burst frequency profiles by considering phaseplane trajectories. In their analysis of downstream models they first plotted the average values of eye velocity versus firing rate for saccades of a given size and direction and noted that the phase plane plot was an open loop. They found that by subtracting the appropriate eye acceleration sensitivity from the BN discharge, the phase plane plot collapsed to a line. Consequently, they concluded that BN activity was related to both eye velocity and acceleration during saccadic eye movements. We believe that the acceleration term in their analysis is the result of two factors: 1) a different frequency content in the neural activity and eye movement signals and 2) an inappropriate choice of lead time.
Our approach differed from that of Van Gisbergen et al. (1981). We used a spike density frequency plot in which the spike train was convolved with a Gaussian pulse whose 5ms width was chosen so that the neural activity was effectively filtered at a comparable frequency to that used for filtering the eye movement signals (Cullen et al. 1996). (To fit a transfer function in a linear system, it is necessary to have the same frequency content in the input and output.)In the latter paper we demonstrated that if the reciprocal interval was used, the firing frequency profile contained frequencies higher than those of the saccade, and this could result in an overestimation of the acceleration term in the downstream models. Furthermore, the reciprocalinterval method is nonlinear and sensitive to noise particularly at high frequencies (Richmond et al. 1987; Sanderson and Kobler 1976). In Van Gisbergen et al. (1981) the BN firing was calculated by using reciprocal intervals, but they averaged the firing rate of BNs during several amplitude matched saccades to reduce the highfrequency noise. It is possible that the significant acceleration sensitivity reported in their analysis could still reflect a difference in the frequency composition of the input eye velocity and output unit activity signals. Moreover, as discussed in Cullen et al. (1996), the choice of lead time strongly influences the estimate of the eye acceleration term. Whereas in our analysis the lead time was determined through optimization, Van Gisbergen and colleagues (1981) used an algorithm in which the lead time was simply defined as the time between the onset of the burst (when, arbitrarily, the discharge rate exceeded 100 spikes/s) and the saccade. As noted in Cullen et al. (1996), we could effectively produce a significant eye acceleration coefficient by shifting, by a short time interval, the unit discharge relative to the corresponding eye movement. A final point worthy of mention is that Van Gisbergen and colleagues (1981) analyzed the difference between the onD and offD responses of BNs, whereas we have chosen to analyze IBN onD and offD responses separately. However, we do not believe that this difference would have greatly effected our characterization of IBNs, because approximately onehalf of our neurons generated negligible offD responses, whereas others were characterized equally well by onD and offD models such that their overall response could be obtained by linearly adding model terms.
Comparison between SLIBNs and LLIBNs
We have considered the discharge of our IBNs in the context of the local feedback model of saccade control (Fig. 2 B) whereby the burst generator is driven by a signal representing motor error [to be analyzed in companion paper III (Cullen and Guitton 1997b)] and in turn drives MNs that move the eye against the resistive forces of the eye plant. In doing so we have explicitly assumed that IBNs and their excitatory counterpart EBNs are similar in their discharge characteristics. This is a common assumption (Fuchs et al. 1985; Hepp and Henn 1983; Luschei and Fuchs 1972; Scudder 1988) and is based on observations that 1) both EBN and IBN discharges encode the amplitude and velocity of saccades (Kaneko et al. 1981; Scudder et al. 1988; Strassman et al. 1986a,b; Van Gisbergen et al. 1981; Yoshida et al. 1982) and 2) EBNs project to the IBN area (Sasaki and Shimazu 1981; Strassman et al. 1986a). An analysis of EBNs, similar to what we have done here, is essential to test this assumption.
An important result of our analysis is that there was no significant difference between SLIBNs and LLIBNs in the VAF by each of the models we tested; i.e., the dynamic signals encoded by SLIBNs and LLIBNs were similar in nature. Hence, although we found a significantly shorter lead time (by 3.7 ms) for SLIBNs, this is not forceful evidence in favor of these cells being closer in their functional connectivity to MNs as is frequently assumed (for review see Fuchs et al. 1985).
Scudder et al. (1988) suggested that 1) LLIBNs might project to omnipause neurons (OPNs) so as to generate the gradual decline in discharge of OPNs before saccades, whereas it is the SLIBNs that provide the pause in antagonist MNs, and 2) LLIBNs may not be inhibitory at all but rather longlead EBNs (LLEBNs) that project to SLIBNs in their vicinity thereby generating the burst in these latter cells. (Our 3.7ms lead time of LLIBNs over SLIBNs is compatible with this scheme.) Our results showing a similarity between SLIBNs and LLIBNs in terms of their dynamic predictions of eye trajectory (Tables 2 and 3) as well as a continuity between the SLIBN and LLIBN populations in terms of their dynamic lead time(s) argue against this categorization of IBNs into different functional populations. Furthermore, there is no experimental evidence showing IBN projections into their own area; Strassman et al. (1986b) have shown, using intracellular staining techniques, that at least some LLIBNs project directly to the motor nucleus. However, the LLIBN projection to contralateral abducens MNs may be functionally weak (Scudder 1988) because LLIBN dynamic lead time begins l5.5 ms before a saccade (i.e., ∼6.5 ms before the pause in the antagonist MNs) and this is very long for a monosynaptic connection. As Scudder et al. (1988) also pointed out, there is no evidence of a decline in MN tonic discharge beginning during the LLIBN prelude. Clearly, more work is needed to determine the role of LLIBNs and LLEBNs in oculomotor control and to distinguish these cells from other longlead BNs (LLBNs) that are implicated in neck motor control (Grantyn and Berthoz 1985).
Conclusions: implications of IBN headfixed model fits
Previous analyses of MN discharges during smooth sinusoidal eye movements have demonstrated that they have a resting discharge and “steady state” eye velocity and acceleration sensitivities (Fuchs et al. 1988; Keller 1973; Van Gisbergen et al. 1981). In addition, these neurons demonstrate a postsaccadic slide in discharge (required by the pole term in Eq. 1 ) (Fuchs et al. 1988; Goldstein and Robinson 1984; Stahl and Simpson 1995). However, the burst discharge of MNs during saccades was not carefully examined using techniques comparable to those that we used here. A thorough characterization of the dynamics of MN activity during saccades would provide an appropriate dynamic model of their discharge, which is important to evaluate the relative contributions of the many premotor inputs to MNs during saccades.
To calculate the transfer function between a single BN (BN_{i}) and the eye movement output, we have pooled the transfer functions together. What is the significance of pooling model parameters and is it a legitimate evaluation of the average BN signal? Van Gisbergen and Van Opstal (1989) have considered this problem with regard to the relationship between ocular MNs and plant time constants. To discuss this point they propose the thought experiment wherein an eye muscle motor unit is paralyzed while its MN is being recorded. The eye movement will surely not be altered by this very minor intervention and thus, they argue, the considerable differences in the time constant values (b _{1}/b _{0} in Eq. 2 ) among individual MNs reflect more differences in inputs to, and intrinsic properties of, MNs rather than the existence of different time constants in the muscle motor units. It follows then that the average b _{1}/b _{0} does not necessarily give the plant time constant.
The same thought experiment can be considered for BNs. Suppose we were to block the axon being studied. The eye movement would undoubtedly remain virtually unchanged because of the large number of BNs acting together. Hence, whereas the transfer function between the recorded cell and the output would remain the same, the neuron is no longer affecting the output. This is because the cell's transfer function is accounting more for the simultaneous input of many BNs onto a MN, and in turn many MNs driving a movement, rather than for its own effect on the movement. It follows that any transfer function is a sample within the IBN population variance itself. However, contrary to MNs, we have found only a small variability in the parameter and VAF values for each IBN and the distribution of parameter values is close to normal. It therefore appears legitimate in our case to take mean values of each parameter and give a mean transfer function for our population.
Figure 18 A reviews inputs to the ABD on the right and left (ABD_{R} and ABD_{L}, respectively) that are implicated in the generation of rightward saccades. Although many cell types have been described in the VN and prepositushypoglossi nucleus (PH), only a few types have been identified as projecting to ABD. The major excitatory input to ABD_{R} from VN_{L} arises from type I positionvestibularpause cells (IPVP_{L}) (cat: Escudero et al. 1992; monkey: Cullen 1991; McCrea et al. 1987; Scudder and Fuchs 1992). During saccades to the right (their onD) most of these cells have tonic activity, some burst, and only a minority pause (rhesus monkey: Scudder and Fuchs 1992). A minor projection arises contralaterally from burst position (_{R}BP_{L}) and position (_{R}P_{L}) cells that are also active during right saccades as indicated by the first subscript. Inhibitory projections from VN to ABD_{R} arise in VN_{R}. However, IPVP_{R}, _{L}BP_{R}, and _{L}P_{R} have their onD to the left and they pause for saccades to the right (– – –). Another cell type in the VN_{R}, the eyehead velocity cells (_{L}EH_{R}), appear to be inhibitory to ABD_{R} (Scudder and Fuchs 1992) but these also pause for rightward saccades.
With regards to the PH, the onD of its neurons is to the ipsilateral side. Excitatory inputs to the ABD_{R} come from cells whose discharge is proportional to eye position to the right (_{R}P_{R}) (Escudero et al. 1992). _{R}BP_{R} in the right VN–PH appear to also send excitatory projections to the ABD_{R} (McCrea et al. 1987). Inhibitory input comes from _{L}P_{L} in PH_{L}, but these neurons pause for rightward saccades.
Let us now consider the projection of EBNs and IBNs to ABD and to the VN–PH, the latter being of considerable interest because of the key role VN–PH plays in the generation of the tonic signal needed to hold the eye at a fixed position in the orbit. As we saw above, for a rightward saccade PH_{R} increases its tonic activity, which is projected to MNs, and PH_{L} decreases its activity. The projection into the PH_{L} of right IBNs (IBN_{R}) is thought responsible for this decrease in tonic activity. Similarly, the projection into the PH_{R} of right EBNs (EBN_{R}) increases the tonic activity. In the simple model of MN discharge given by Eq. 2 , it is assumed that the tonic signal is obtained by integrating the BN eye velocity signal (reviewed in Fuchs et al. 1985; Robinson 1975). VN–PH classically has been thought of as this neural integrator (Fig. 2 B) (see Cannon and Robinson 1987; Cheron et al. 1986a,b; Cheron and Godaux 1987; Escudero et al. 1992; LopezBarneo et al. 1982; McFarland and Fuchs 1992; Mettens et al. 1994). Our conclusion that IBN firing frequency is not simply proportional to E˙ complicates the interpretation of how the integration process is performed.
The net signal of BNs onto ABD and VN–PH is the combination of the onD and offD discharges of EBNs and IBNs. For a saccade to the right, the right IBN population produces the signal
In conclusion, our analysis using rigorous system identification algorithms reveals that understanding the signal processing in even a “simple” structure like the ABD of the headfixed monkey is a very complex and laborintensive process. We cannot help but be struck by the awesome task of resolving the much more complex circuits involved in mediating higher brain functions.
Acknowledgments
We are grateful to C. G. Rey for useful discussions and contributions to the development of the analysis procedure and software and to H. L. Galiana for helpful discussions. We thank Dr. J.A.M. Van Gisbergen for a number of useful suggestions that helped to improve this manuscript greatly. We also thank Dr. W. M. King for critically reading this manuscript.
This study was supported by the Medical Research Council of Canada, the National Institutes of Health, and the Human Frontiers Science Organization.
Footnotes

Address for reprint requests: K. E. Cullen, Aerospace Medical Research Unit, 3655 Drummond St., Montreal, Quebec H3G 1Y6, Canada.
 Copyright © 1997 the American Physiological Society