## Abstract

We compared the spike activity of individual neurons in the*Aplysia* abdominal ganglion with the movement of the gill during the gill-withdrawal reflex. We discriminated four populations that collectively encompass approximately half of the active neurons in the ganglion: (1) second-order sensory neurons that respond to the onset and offset of stimulation of the gill and are active before the movement starts; (2) neurons whose activity is correlated with the position of the gill and typically have a tonic output during gill withdrawal; (3) neurons whose activity is correlated with the velocity of the movement and typically fire in a phasic manner; and (4) neurons whose activity is correlated with both position and velocity. A reliable prediction of the position of the gill is achieved only with the combined output of 15–20 neurons, whereas a reliable prediction of the velocity depends on the combined output of 40 or more cells.

Neural systems can be viewed as complex input–output devices. The input can be an external stimulus, and the output is the generated behavior. From the point of view of neuronal computation, two essential and related questions emerge. First, how is the sensory input and intended motor output represented? Second, how distributed is the processing in the neural system?

To help answer the above questions, we had monitored the activity of a large fraction of the population of the neurons that generated a behavioral response. Here we determined how well the behavior could be both fit and predicted with weighted summations of the spike activity of the individual neurons. In particular, we analyzed the spike activity of 149 neurons in the abdominal ganglion of *Aplysia*and the simultaneously recorded gill movements that occurred after stimulation of the siphon with a light mechanical touch (Wu et al., 1994). The spike data were obtained from voltage-sensitive dye recordings that allowed simultaneous monitoring of the activity of a substantial fraction of the active population of neurons. These neurons are thought to be mainly interneurons and motor neurons because the primary sensory neurons for light touch are probably in the periphery and the sensory neurons in the ganglion have a very modest response to the mechanical stimulus that was used in these experiments (Hickie et al., 1997).

It had been suggested that there is a simple functional architecture for the gill-withdrawal reflex. Specifically, Byrne et al. (1978)estimated that the feedforward circuit formed by monosynaptic connections between eight LE sensory neurons and six gill motor neurons can account for 60% of the motor neuron postsynaptic potential. However, more recent evidence shows that the contribution of the monosynaptic LE sensory component to the movement is actually an order of magnitude smaller, ∼5 versus 60%, and that contributions from other sensory neurons as well as interneurons are important (Hawkins et al., 1981; Frost et al., 1988; Cohen et al., 1991; Trudeau and Castellucci, 1992; Hickie et al., 1997; Walters and Cohen, 1997). In addition, voltage-sensitive dye measurements (Zecevic et al., 1989;Nakashima et al., 1992; Tsau et al., 1994; Hopp et al., 1996) suggest that ∼300 of the ∼1000 neurons in the ganglion are activated during gill-withdrawal reflex; this allows the possibility that the coding of gill movement may be distributed among many neurons.

Here we analyze the previously recorded spike data, and ask the following. (1) Are different aspects of the gill-withdrawal movement, position as opposed to velocity, independently coded by neurons in the abdominal ganglion? (2) If different aspects are independently coded, does this involve overlapping or separate pools of neurons? (3) How large a population of neurons in the abdominal ganglion is required to reliably predict the movement of the gill based on the weighted output of neuronal activity?

## MATERIALS AND METHODS

The experimental data for spike activity that we analyzed were obtained earlier by Wu et al. (1994). The experiments were performed on an isolated siphon preparation developed by Kupfermann et al. (1971). The data consist of recordings made during seven separate light mechanical touches, each consisting of a force of 10 mN that was applied to the siphon for 400 msec. The interstimulus interval was 15 min and was chosen to minimize habituation. The neuronal activity in the ganglion was determined by measuring the light transmitted through a ganglion stained with the voltage-sensitive oxonol dye JPW1131 (ne RH155) (Grinvald et al., 1980). Because of the limited signal-to-noise ratio in the measurements, it was estimated that only half of the active neurons in the ganglion were detected (Wu et al., 1994). The behavioral data consists of a time series that represents the area of the gill. The *top* of Figure 1shows a portion of the spike recordings, before as well as after the siphon stimulation (Fig. 1, *S*), and the *bottom*shows the recordings of the gill area; these areas were remeasured for this study to improve the temporal accuracy. The data from trial 8 in the data of Wu et al. (1994) was omitted because there were fewer spikes during that trial and subsequent analysis of that trial was numerically unstable.

*Analysis.* We have both the time series representing the gill movement and the simultaneously recorded spike data from ∼150 neurons in seven consecutive trials. We attempted to reconstruct the gill movement from the spike times of the individual neurons in two different ways. First, we fit the behavior from each trial from the spikes times of that trial in terms of a set of weights that parameterized the fit. This allowed us to examine the form and consistency of the weights derived from independent trials. Specifically, we attempted to fit the gill movement in terms of the optimal linear combination of the output of each neuron. The weights used for each neuron, between −1.0 and +1.0, are the free parameters of the fit; these weights are adjusted by an algorithm (see below) that minimizes the difference between the recorded movement and the fitted pattern.

As a second measure of the relationship between the spike times and behavior, we predicted the behavior of one trial using a single set of synaptic weights fitted as an average over the other six trials. If the neuron weights are consistent from trial to trial, then the behavioral curve of the one remaining trial should be reasonably well predicted.

We made several simplifying assumptions that were common to both analyses. First, we assumed that each spike in a train would make an equal contribution to the behavior, and thus, we did not consider synaptic facilitation or depression. Second, we assumed that there was a linear relationship between spikes in neurons and the position of the gill and the velocity of the gill movement, and thus, we did not consider possible interactions between neurons. Last, because the spikes are discrete events and the behavioral response is a continuous function of time, we represented each spike by a gaussian curve with an SD of 70 msec. After this transformation, the output activity of each neuron is represented by a sequence of gaussian time courses.

*Algorithms.* We define the predicted behavioral response, either the position or the velocity of the gill, to be a linear combination of neural activity during that trial. Formally, this is given by:
Equation 1where *k* is the trial number, *N* is the number of neurons considered for the fit (*n* = 133 for the cells included) and the *w _{i}
* are the weights that describe the relative influence of individual neurons on the behavioral pattern. Those weights are the free parameters of the fit.

*Spike train representation.* The time series*S ^{k}
_{i} *represent the temporal activity of a single neuron in terms of its spike times. Ideally, because the spikes are discrete events and the behavioral response is a continuous function of time, one would smooth each spike by convoluting it with a postsynaptic response function. However, as a simplifying procedure that preserves the essence of our analysis, every spike is converted to a normalized gaussian with the width of ς that is centered at the time of the spike, i.e.,
Equation 2where

*t*denotes time,

*s*denotes the spike train of

_{ik}*i*-th neuron in

*k*-th trial, and

*t*is the time of a single spike in that spike train. The gaussian smoothing may be viewed as assigning a probability distribution that a given neuron fired at a particular time, under the assumption that this probability is the highest at the time the neuron

_{ik}*de facto*fired and is symmetrical. The SD was set to ς = 70 msec for both velocity and position fits. Qualitatively similar results were obtained with values for ς between 30 and 200 msec.

Because the changes in the position component are much slower then 70 msec used for the fits, the fits of this component were additionally smoothed with a sliding window of 1000 msec width. The smoothed fits (see Fig. 4) were used in the calculation of the error of fit.

*Individual trial fits.* The aim of the method is to optimize the weights of the neurons in such a way that the behavior reconstructed from the spike series will best match the position or velocity during the withdrawal reflex. This is done through the standard procedure (Korn and Korn, 1968) of minimizing a quadratic error function, *E _{k}
*,
Equation 3where

*B*denotes the recorded behavior, either position or velocity of the gill,

^{obs}_{k}*B*is the fit that is being optimized, and

^{pre}_{k}*t*= 0 denotes the onset of stimulation whereas

*T*= 10 sec is the time at the end of the record. The function

*E*measures the difference between predicted and observed behavior. Our aim is to find the set of weights that minimizes the error function on a given trial. Those weights are be then used to reconstruct the movement.

_{k}The values of the weights that minimize the error function are obtained directly from:
Equation 4which has to be fulfilled to find the minimum of function*E _{k}
* by varying the weights

*w*; the arrow in Equation 4 denotes a vector of the weights of all the neurons. We noticed that the solution to this equation can be written in the analytical form: Equation 5where Φ is a matrix defined by: Equation 6The function Φ represents the overlap of activity of the spikes from the

*i*-th and

*j*-th neurons. It estimates the relative importance of the given neuron. If there are many neurons that fire in a similar pattern, the relative importance of the single cell from that group is smaller; on the other hand, if there is a single cell with a given pattern of activity, its relative importance to the fit may be higher. Thus, if two neurons spike at the same time, the exponent achieves its maximal value. If the spike times do not match, the exponent will tend to 0. The weights are proportional to the inverse of this matrix.

The matrix Φ is sparse because of the relatively low spike rate of most cells. In particular, only 24 of the 133 neurons had more than three spikes in all seven trials, and only 13 of the 133 had more then five. This sparseness may cause Φ to be singular and is expected to lead to numerical instabilities in the algorithm during the matrix inversion. To calculate the inverse of matrix Φ under such conditions, we performed a spectral decomposition (Golub and Van Loan, 1996) and used the leading eigenvectors from that decomposition to estimate the inverse of Φ.

The function Ψ is defined as: Equation 7and measures the correlation of the recorded position or velocity component of the movement with the spike trains of a given neuron. The position and velocity were rescaled in such a way that the values of the position and velocity before the stimulus were considered as the baseline and were set to 0. All trials and fits were normalized, so that the position had maximum value of 1. Note that the weights were obtained through numerical solution of Equation 5 and not through a numerical search for minima of the error function surface extended in the space of all possible weights. This avoided the problem of converging on local minima. The calculations were written and performed in Interactive Data Language (Research Systems Inc, Boulder, CO).

*Optimal weights for novel trial predictions.* To predict the movement in one trial based on the measurements for the other six trials, either the position or velocity, we performed an analogous procedure as described above. However, the error function is now defined so that the sum of errors of the six other trials is, i.e.,
Equation 8The weights that result from this minimization are then used to reconstruct the behavioral curve, position or velocity, of the remaining trial. This procedure is repeated for every permutation of trials. As before, Equation 5 is an analytical solution to the minimization problem with the vector Ψ and the matrix Φ now given by:
Equation 9
Equation 10

## RESULTS

The neuronal data (Fig. 1, *top*) suggests a possible dichotomy between two populations of neurons. As shown at the*top* of the raster diagram, there is a small population of cells that respond to the stimulus by firing tonically, with a slow modulation. On the other hand, there is a large population of neurons that fire infrequently but to some extent synchronously. This leads to the speculation that the two groups of cells might correlate with aspects of the gill movement that occur on different time scales. In particular, the recorded movement of the gill (Fig. 1,*bottom*) can be characterized by two dynamic components: (1) slow contractions and relaxations that can be described on a time scale of 500–1000 msec; and (2) fast contractions that always occurred at the beginning of the gill withdrawal and sometimes occurred later in the response. The relevant time scale for the fast contractions is ∼100 msec.

To quantify the two dynamic components, we decomposed the movement data, as illustrated for trials 1 and 5 (Fig.2*A*). (In this and subsequent figures, the onset of the mechanical stimulus is the 0 time.) Gill position is found by filtering the recorded movement with a low-pass filter (cutoff frequency of 0.6 Hz) to remove the fast contractions from the recorded movement (Fig. 2*B*); the peaks of gill contraction were not shifted compared with the original traces. The fast contractions are equated with the velocity of the gill movement, given by the first derivative of the movement (Fig.2*C*). In addition, Figure 2*D* shows the two cumulative histograms of the spike activity of all the neurons from Figure 1 for the two trials. Note that there are peaks in the histograms that are coincident with peaks in the velocity traces (Fig.2*C*). This coincidence, combined with the results in Figure1, suggests that there is a subpopulation of cells whose activity correlates with velocity. This hypothesis is quantitatively analyzed below.

### Second-order sensory neurons

A subpopulation of “sensory” neurons was defined by identifying neurons that spiked during the period of stimulus presentation (400 msec) but before the movement started. To avoid including neurons that fired during the stimulation period by chance, we added the constraint that the neuron had to have at least five spikes in the seven trials during this period. This group consisted of 16 neurons. Figure 3*B* shows the cumulative histograms of the activity of these 16 cells for two trials, as well as the velocity traces for those two trials. In both trials, the first peak in the histogram comes before the onset of the movement. There is a second peak that, in some instances, is synchronized with a velocity peak. However, previous results from habituated preparations, in which there is no gill movement, strongly suggests that the second burst is the response to the offset of the mechanical stimulus (Falk et al., 1993). The size, location in the ganglion, and sensitivity to altered Ca^{2+}concentrations of this “sensory” class makes it likely that they are second-order sensory neurons (Hickie et al., 1997). We removed the group of 16 putative second-order sensory neurons from computations shown in this paper.

### Fast and slow gill movements correlate with activity of different neural populations

We fit the gill position and velocity with the neural activity in each of the seven trials separately. Examples of the gill position and velocity fits from two trials are presented in Figure4*A*; the gill position and velocity are shown in *gray*, and the fits are *thin dark lines*. The entire neural population, except the second-order sensory neurons, was used to calculate both fits. Comparison of the behavior and the fits shows that both velocity and position are well reproduced by the fitting algorithm (Eqs. 5-7), indicating that both position and velocity have their cellular correlates.

We quantify the goodness of the fit in terms of a “percentage error,” denoted *PE _{k}
*. This is a normalized measure, defined as the mean-square distance between the position (or velocity) traces and their respective fits (Eq. 3) divided by the mean-square amplitude of the position (or velocity), times 100, i.e.,
Equation 11where

*B(t)*refers to the position or the velocity component of the behavior (Fig.2

*B*,

*C*). The percentage error is

*PE*= 0% if the fit is perfect, and it tends toward 100% as the fit worsens. The percentage error for the two position fits in Figure 4 are

_{k}*PE*

_{1}= 4% and

*PE*

_{5}= 6%. The percentage error for the two velocity fits in Figure 4 are

*PE*

_{1}= 13% and

*PE*

_{5}

*=*28%. The trial average percentage error (Eq. 11 with

*PE*averaged over

_{k}*k*) was 5 ± 1% (mean ± SD of the mean) for all seven position fits and 25 ± 6% for the velocity fits.

To determine whether the different aspects of the gill withdrawal are correlated with separate or overlapping populations of neurons, we examined the distribution of the weights for the two fits. Figure5 shows the weights of the 32 neurons assigned the largest weights (‖*w*‖ > 0.25 for the fit to either gill position or velocity); the mean value of each weight is an average over all trials of the separate fits to position and velocity, and the bars denote the SD. We observe that approximately half of the cells that correlate with position, i.e., which have a weight whose magnitude is large when fitted to the position traces, do not correlate with velocity (Fig. 5, *red symbols*). Furthermore, there is a group of neurons that correlates with velocity but not with position (Fig. 5, *blue symbols*). Thus, the major fraction of neurons exclusively correlates with only one component of the motion. A substantial group of neurons, however, also correlates with both position and velocity. This latter population is distinguished by containing both positive and negative weights for the velocity component.

It is instructive to examine the firing patterns of the populations of neurons as a means to understand their differences in output in relation to the movement of the gill. Figure 4*B* shows the activity of the six neurons with the largest positive weights for position and velocity for trials 1 and 5. Some of the neurons whose activity is positively correlated with position have many spikes, and their cumulative activity tends to be distributed over the whole interval of the response. On the other hand, the positive weighted neurons whose activity is correlated with the velocity are phasic and fire much less frequently albeit somewhat synchronously.

Although neurons that significantly contribute to the position do so almost exclusively through positive weights (Fig. 5), the velocity receives a significant contribution from negatively weighted neurons (Fig. 5). For example, in both trials 1 and 5, one of the tonically firing neurons that has large positive weight for position also has a large negative weight for velocity (Fig. 4*C*, ★ ★). This neuron was silent when the velocity was large; an increase in velocity may correspond to a release from inhibition from this neuron.

### Position and velocity correlates are different

To quantify the above differences between the activity of the neurons driving position and velocity, we measured the entropy of the cumulative spike histograms of the 20 largest positively weighted cells in the two fits. The entropy quantifies the degree of randomness of the distribution of neuronal spiking over time and is defined as:
Equation 12where *p _{I}
* is the probability of the event at the given time and

*I*labels the time bin, whose width is taken to be the same as that for the gaussian broadening spike, i.e., ς = 70 msec (Eq. 2). The average entropy for the group of neurons used to fit position was

*S*= 2.4 ± 0.1 bits (mean ± SEM), whereas for the group of neurons used to fit velocity the value was

*S*= 2.0 ± 0.1 bits. This difference, 0.4 bits, is statistically significant (

*p*= 0.001) and is consistent with the observation that most of the neurons that generate the velocity component fire phasically and to some extent synchronously, whereas many of the neurons whose activity correlate with the position fire spikes tonically and in an incoherent manner.

### Trial-averaged optimal weights and prediction of the behavior

To investigate the consistency of the neuronal populations that correlate with position and velocity, we used the weights from six of the trials as a means to reconstruct the remaining trial (Eq. 1). Specifically, we calculated the weights from a minimization procedure that used activity and movement patterns of the six trials simultaneously (Eqs. 8-10). Those weights, in conjunction with spike activity in the seventh trial, were used to construct the behavioral curves, both position and velocity, for that trial. If the activity of the same population consistently correlates with a given aspect of the movement, the predicted curve will resemble the recorded one. The result of the predictions with trial 1 as the test trial and trial 5 as the test trial are presented in Figure 6. For both trials and for both position and velocity, the predictions fit the recorded behaviors well, although not as well as the fits of individual trials (Fig. 4*A*). The *PE* values for the two predictions in Figure 6 were*PE*_{1} = 5% and*PE*_{5} = 8%, respectively, for position and *PE*_{1} = 44% and*PE*_{5} = 40%, respectively, for velocity. The match between the measured and predicted movements of the gill was calculated for all seven possible permutations. Critically, the prediction compared favorably with the measured behavior in all trials, although the behaviors can be quite different from trial to trial (Fig. 1). The average percentage error (Eq. 11 with*PE _{k}
* averaged over

*k*) for all seven possible predictions was 9 ± 4% for position and 43 ± 10% for velocity, also somewhat larger than the percentage error for the individual fits.

### Population size

As a means to assess the size of the subpopulations whose activity correlates with position and velocity, we restricted the number of neurons that were used for the predictive fits and examined the effects of these restrictions on the accuracy of the predictions. The number of cells was limited in either of two ways: (1) we used only a limited number of the “best” neurons, i.e., those with the largest absolute weights; or (2) we omitted a number of best neurons. For both position and velocity, we computed the percentage error after limiting the neuron populations to the 2, 5, 10, 15, 20, 30, 40, 50, and 60 best, and after omitting the 10, 20, 30, 40, 50, and 60 best. The *left side* of Figure 7 presents the percentage error (Eq. 11 with *PE _{k}
* averaged over

*k*) for the cases when the activity of a limited number of best neurons is included in the predictions; the

*right side*has the percentage error for the cases when the best cells are omitted from the predictive fit.

The results indicate that the number of position-correlated cells required to achieve optimal prediction is ∼20. The percentage error for prediction of position worsens by one SD when we limit the number of cells used for the prediction below ∼15 (Fig. 7, *left*). Furthermore, when the most strongly weighted cells are omitted from the predictions (Fig. 7, *right*), the percentage error for the position increases quickly, even if only 10 or 20 such cells are omitted. In contrast to the case for position, the number of velocity-correlated cells appears to be much larger. The prediction improves monotonically when up to 60 neurons are included in the fit and worsens monotonically when a increased number of such strongly weighted cells is omitted.

As a control to determine the extent to which the fitting algorithm can fit spurious data, we repeated the procedure after shifting the behavior by up to 4.0 sec. The asymptotic value of the errors of the predictions for the shifted data (using all of the neurons) are shown as the *horizontal dashed black* and *gray lines* on Figure 7. For the shifted data, the errors are large for both position and velocity.

## DISCUSSION

We used optical recordings from the *Aplysia* abdominal ganglion (Fig. 1) and a numerical algorithm (Eqs. 1-10) to discriminate functionally different neuronal populations involved in the gill-withdrawal movement. Our results indicate that the gill-withdrawal reflex consists of two elements (Fig. 2): (1) relatively slow contractions and relaxations, i.e., position; and (2) relatively fast contractions, i.e., velocity. These two elements, which form one continuous behavioral pattern, are correlated with two partially separate neural populations (Fig. 5) that appear to use two different coding schemes. Position correlates, in large part, with a relatively small population of tonically active neurons whose activity is slowly modulated throughout the response (Fig. 4). These cells may code for position via excitation. Velocity correlates primarily with a relatively large population of infrequently firing cells that generate narrow bursts of activity (Fig. 4). These cells may code for velocity via excitation. In addition, velocity is also correlated with a few tonically firing cells that cease firing for brief epochs during velocity peaks (Figs. 4, 5); in this case, an increase in velocity may result from a release from inhibition.

In the group of 133 neurons that we included in the fits and predictions, the output from ∼15–20 cells are needed to reliably predict the position of the gill (Fig. 7). The output activity of a larger population, at least 40 neurons, is needed to predict velocity (Fig. 7). Earlier estimates of the completeness of the optical recordings (Wu et al., 1994) suggested that only half of the active neurons were detected. Assuming that the other half of the neurons have the same distribution of properties as the detected half, then the best estimate for the number of neurons correlated with position is at least 30, and the number correlated with velocity is at least 80. We emphasize that our results do not establish a causative relationship between the activity of recorded neurons and the movements itself. For example, we do not know whether the neurons we found are interneurons or motor neurons.

The firing patterns of the previously identified motor neurons (Kupfermann and Kandel, 1969; Peretz, 1969; Carew et al., 1974;Kupfermann et al., 1974; Koester and Kandel, 1977) suggest that they are tonically firing neurons of the type that received large positive weights in the position fit. On the other hand, the number of neurons that control the velocity component is large, and many of the neurons fire very sparsely. Thus, their individual contribution to gill contraction would have been difficult to detect in microelectrode studies. Their number and role becomes apparent when they can be monitored as a large population of cells.

Our results indicate that different components of the gill-withdrawal movement correlate with neuronal activity patterns using two different coding schemes: average rate coding of frequently firing cells for position, and coding primarily through coincident activity of a large population of infrequently firing cells for velocity. Results of experiments on mammalian motor systems at the level of either brainstem nuclei (Precht 1979) or neocortex (Georgopoulos et al., 1986, 1999;Schwartz, 1992; Moran and Schwartz, 1999) suggest that few cells are specific to position or velocity, but instead each cell codes for both parameters. Our results, however, indicate that, in *Aplysia*, the neuronal populations that are correlated with velocity and position are partially separate. The difference between our imaging-based results for *Aplysia* and the electrode-based results for cortex could be attributed to sampling bias, e.g., the infrequent spiking of most neurons in the velocity-correlated population could allow them to be overlooked during electrode recordings.

Our results show that the activity of large populations of neurons is needed to fit or predict both position and velocity. Additionally, we found (Fig. 7) that limiting the number of cells included in the prediction reduced the quality of the prediction, which indicates that many cells are needed to reliably recreate the behavioral curves. This implies a distributed coding of the gill-withdrawal reflex, at least in the sense that the output from one or few neurons has insufficient reliability to accurately code the movement of the gill. In sensory areas of mammalian cortex, reliable decoding of sensory input typically requires averaging over the activity of many neurons (Seung and Sompolinsky, 1993), although examples of reliable coding by single cortical cells have been shown (Newsome et al., 1989; Fee et al., 1997). In further concurrence with the results for gill withdrawal in*Aplysia*, only distributed processing has been identified in the context of motor control by mammalian cortex (Georgopoulos, 1995).

## Footnotes

This project was initiated as a part of the National Institute of Mental Health sponsored course “Methods in Computational Neuroscience” at the Marine Biological Laboratory. The work was also sponsored in part by National Institutes of Health Grant NS08437, a Brown-Coxe fellowship from the Yale University School of Medicine, and the Grass Foundation. We thank William Bialek, David Hansel, Bill Ross, Haim Sompolinsky, Terry Walters, and Dejan Zecevic for discussions. We also thank the reviewers for several helpful suggestions.

Correspondence should be addressed to Michal Zochowski, Department of Molecular and Cellular Physiology, Yale University School of Medicine, 333 Cedar Street, New Haven, CT 06520. E-mail: mrz{at}fred.med.yale.edu.