The Journal of Neurophysiology Vol. 79 No. 2 February 1998, pp. 1017-1044
Copyright ©1998 by the American Physiological Society

Interpreting Neuronal Population Activity by Reconstruction: Unified Framework With Application to Hippocampal Place Cells

Kechen Zhang¹, Iris Ginzburg¹, Bruce L. McNaughton², and Terrence J. Sejnowski^{1, 3}

¹ Computational Neurobiology Laboratory, Howard Hughes Medical Institute, The Salk Institute for Biological Studies, La Jolla, 92037; ² Division of Neural Systems, Memory, and Aging and Department of Psychology, Arizona Research Laboratories, University of Arizona, Tucson, Arizona 85724; and ³ Department of Biology, University of California, San Diego, La Jolla, California 92093

	ABSTRACT

Abstract Introduction References

Zhang, Kechen, Iris Ginzburg, Bruce L. McNaughton, and Terrence J. Sejnowski. Interpreting neuronal population activity by reconstruction: unified framework with application to hippocampal place cells. J. Neurophysiol. 79: 1017-1044, 1998. Physical variables such as the orientation of a line in the visual field or the location of the body in space are coded as activity levels in populations of neurons. Reconstruction or decoding is an inverse problem in which the physical variables are estimated from observed neural activity. Reconstruction is useful first in quantifying how much information about the physical variables is present in the population and, second, in providing insight into how the brain might use distributed representations in solving related computational problems such as visual object recognition and spatial navigation. Two classes of reconstruction methods, namely, probabilistic or Bayesian methods and basis function methods, are discussed. They include important existing methods as special cases, such as population vector coding, optimal linear estimation, and template matching. As a representative example for the reconstruction problem, different methods were applied to multi-electrode spike train data from hippocampal place cells in freely moving rats. The reconstruction accuracy of the trajectories of the rats was compared for the different methods. Bayesian methods were especially accurate when a continuity constraint was enforced, and the best errors were within a factor of two of the information-theoretic limit on how accurate any reconstruction can be and were comparable with the intrinsic experimental errors in position tracking. In addition, the reconstruction analysis uncovered some interesting aspects of place cell activity, such as the tendency for erratic jumps of the reconstructed trajectory when the animal stopped running. In general, the theoretical values of the minimal achievable reconstruction errors quantify how accurately a physical variable is encoded in the neuronal population in the sense of mean square error, regardless of the method used for reading out the information. One related result is that the theoretical accuracy is independent of the width of the Gaussian tuning function only in two dimensions. Finally, all the reconstruction methods considered in this paper can be implemented by a unified neural network architecture, which the brain feasibly could use to solve related problems.

	INTRODUCTION

Abstract Introduction References

Reconstruction is a useful strategy for analyzing data recorded from populations of neurons, in which external physical variables such as the orientation of a light bar on a screen, the direction of hand movement in space, or the position of a freely moving animal in space are estimated from brain activity. Reconstruction is sometimes called decoding, where "coding" describes the tuning of neuronal responses to the variables of interest. A wide range of neural coding and decoding problems have been studied previously (Abbott 1994; Bialek et al. 1991; Optican and Richmond 1987; Salinas and Abbott 1994; Seung and Sompolinsky 1993; Snippe 1996; Zemel et al. 1997; Zohary et al. 1994).

Two main goals for reconstruction are approached in this paper. The first goal is technical and is exemplified by the population vector method applied to motor cortical activities during various reaching tasks (Georgopoulos et al. 1986, 1989; Schwartz 1994) and the template matching method applied to disparity selective cells in the visual cortex (Lehky and Sejnowski 1990) and hippocampal place cells during rapid learning of place fields in a novel environment (Wilson and McNaughton 1993). In these examples, reconstruction extracts information from noisy neuronal population activity and transforms it to a more explicit and convenient representation of movement and position. In this paper, various reconstruction methods that are theoretically optimal under different frameworks are considered; the ultimate theoretical limits on the best achievable accuracy for all possible methods also are derived.

Our second goal for reconstruction is biological. Because the brain extracts information distributed among the activity of populations of neurons to solve various computational problems, the question arises as to which algorithms feasibly might be used in the brain. Contrary to previous views, the various reconstruction methods, including the Bayesian methods, can be implemented by a simple feedforward neural network. This demonstrates that biological systems have the resources needed to implement the most efficient reconstruction algorithms, which can reach the Cramér-Rao lower bound under suitable conditions. In addition, our theoretical results on minimal reconstruction error become immediately relevant for understanding the accuracy of population coding in biological systems and its dependence on the tuning parameters of individual cells.

In this paper, hippocampal place cells serve as a primary example of the reconstruction problem and are used for testing different methods. Some interesting biological properties of place cells are revealed by the reconstruction, which illustrates the power of this approach for studying populations of neurons.

	DESCRIPTION OF PLACE CELL FIRING

Representative example of reconstruction

Position reconstruction based on hippocampal place cell activity provides an excellent example that captures all major aspects of reconstruction (Wilson and McNaughton 1993). The task is to infer the position of the rat's head based on the patterns of spiking activity from simultaneously recorded place cells. As illustrated in Fig. 1, a place cell is often quiet and fires maximally only when the animal's head is within certain restricted regions in the environment, called the place fields (McNaughton et al. 1983; Muller et al. 1987; O'Keefe and Dostrovsky 1971; Olton et al. 1978).

View larger version (44K): [in this window] [in a new window]   FIG. 1. Firing rate maps of 25 hippocampal place cells simultaneously recorded in a rat (animal 1) running on elevated track of a figure-8 maze, with clockwise movement in upper loop and counterclockwise movement in lower loop, and slowing down at corners for food. Restricted regions with high firing rates are called place fields. Another 32 simultaneously recorded cells in same animal were active during sleep but virtually silent on this maze and therefore were not shown. Maps were computed from 7 min of continuous data. Reciprocal fields were obtained from firing rate maps by using a pseudoinverse. In each plot, scaling is linear with zero value corresponding to 0 in color map and maximum positive value corresponding to 1. Here reciprocal fields roughly resemble place fields, but have negative values in surrounding regions.

Successful simultaneous recordings from many cells in the hippocampus of a freely behaving rat make it possible to study population activity and reconstruction in single trials. In most other systems, such as recordings from motor cortex during reaching, neurons usually are recorded sequentially during the same task, and the population is examined later under the assumption that the cells would have similar statistics if recorded simultaneously. In single-trial reconstruction, phenomena such as erratic jumps in reconstructed trajectory easily can become apparent (Fig. 9).

View larger version (74K): [in this window] [in a new window]   FIG. 9. Sudden jumps of reconstruction error (here showing Bayesian 1-step method for animal 1) often occurred when instantaneous firing rate averaged over all cells was low, which in turn was correlated with low running speed. In firing rate diagram, dotted curve shows inverses of average firing rates; this serves only as a visual guide to emphasize moments of low firing rates. Running speed and reconstruction error also had correlates with hippocampal local field potentials, because sharp waves and theta rhythm are known to be more prominent during immobility and motion, respectively. Amplitude diagram shows absolute value of local field potential signal of channel 1, which was slightly blurred and shown in arbitrary units. Continuous Fourier power spectra were obtained by sliding a Hanning window with a length 10 times the period of each frequency. Channel 1 was sampled at 1 kHz to capture sharp waves, which can be seen in Fourier power spectra as thin vertical strips ~150-200 Hz. Channel 2 was sampled at 200 Hz to capture theta rhythm, which can be seen in Fourier power spectra as horizontal dark bands ~8 Hz with 2nd harmonics around 16 Hz.

It is important to realize that the general formulation of the reconstruction problem does not rely on any unique properties of place cells. In other words, the same approach used for place cells could be applied to other reconstruction problems without essential modification. The vector x, which is interpreted in these recordings as the position of the animal in the maze, could in general be interpreted as any external physical variables with which neural activity is correlated.

Place cell data

The experimental data analyzed in this paper were obtained by methods that have been described by Wilson and McNaughton (1993), Skaggs, McNaughton, Wilson, and Barnes (1996), and Barnes, Suster, Shen, and McNaughton (1997). The accuracy of spike timing data of neurons from areas CA1 or CA3 of the right hippocampus was 0.1 ms. The head position was tracked by infrared diodes on the head stage that were sampled at 20 Hz and at a resolution of 256 × 256 pixels (~111 × 111 cm). Several channels of hippocampal local field potentials were recorded at 200 or 1,000 Hz.

Data were obtained from two young rats while running continuously in familiar mazes. Animal 1 was running on the elevated tracks of a figure-8 maze (Fig. 1), with clockwise movement in the upper loop and counterclockwise movement in the lower loop. The outer dimensions of the maze were 93.5 × 80.7 cm, and the width of the tracks was 6.4 cm (Fig. 1). Animal 2 was running on the tracks of a rectangular maze, physically identical to the upper loop of the figure-8 maze. Both animals stopped briefly at the corners of the mazes for food reward during the recordings.

In reconstruction, the first 7 min of data from animal 1 were used to compute firing rate maps and visiting probabilities, and the subsequent 7 min of data were used for reconstruction. For animal 2, the sampling and testing intervals were 10 min each. For comparison, we also used the second halves of data for sampling and the first halves for testing. In another analysis, the data were divided into consecutive 3-min time blocks, with only a single block chosen for sampling, leaving all remaining blocks for testing.

Basic description

The goal of this section is to quantify the basic tuning property of a place cell. Reconstruction only needs two distribution functions derived from the data. The first one is the probability P(x) for the animal to visit each spatial position x = (x, y), a vector that denotes the position of the animal's head on the horizontal plane. In our dataset, x was discretized on 256 × 256 grid and sampled at 20 Hz. Let N(x) be the number of times the animal was found at position x during a time period, say, of several minutes; then the probability of finding the animal at x can be estimated as

(1)

This probability is sometimes called spatial occupancy.

The second distribution needed by reconstruction is the average firing rate f(x) of a place cell for each position x. The firing rate distribution f(x) is sometimes called the tuning function or the firing rate map (Fig. 1). Let S(x) be the total number of spikes collected while the animal was at location x, then the average firing rate at x is estimated as

(2)

where t is the time interval of position tracking so that N(x)t is the total time spent at position x. The derived firing rate map f(x) was smoothed by convolving it with a Gaussian kernel (Parzen window) with a standard deviation of 1 cm. In theory, the time average of the firing rate should equal to the spatial average f_mean = f(x)P(x). This is no longer exactly true after the aforementioned smoothing, and therefore can be used as a cross-check.

The firing rate distribution f(x) can be interpreted in terms of conditional probability as follows

(3)

This is equivalent to Eq. 2, as is verified readily by using Eq. 1 and the definition of two additional quantities, namely, the conditional probability

(4)

for an arbitrarily chosen spike to fall into the spatial position x, and the overall average firing rate

(5)

with S(x) being the total number of spikes and N(x)t being the total time of recording. The firing rate distribution f(x) describes the intrinsic tendency for the cell to fire at each position x and is independent of how often the animal visits this position.

Our method of computing the firing rate map is similar to that in Muller, Kubie, and Ranck (1987). An alternative method with adaptive binning has been used by Skaggs, McNaughton, Wilson, and Barnes (1996). Additional consideration of the heading direction and the intrinsic directionality of place cells (McNaughton et al. 1983; Muller et al. 1994) may be formulated by running velocity modulation of firing rate as in the next section.

Firing rate modulation

The description in the preceding section can be broadened to incorporate the modulation of firing rate by additional variables other than the position x, such as running speed, the heading direction, and the theta rhythm. For our dataset, including firing rate modulation only slightly improved reconstruction accuracy.

RUNNING VELOCITY. The firing rate of a place cell typically increases with running speed and also depends on the heading direction when movement is confined to a narrow track (McNaughton et al. 1983). Both the effects of speed and direction can be modeled as multiplicative modulation of firing rate by the running velocity. For our dataset, the animals run mostly in one direction, and we need only to consider the firing rate modulation by speed.

In general, suppose the average firing rate f(x, v) of a cell at position x also depends on the instantaneous running velocity v. Consider multiplicatively separable firing rate modulation of the form

(6)

where functions F(x) and G(v) depend on x and v only, respectively. The simplest case is linear speed modulation

(7)

That is, G(v) = a + b is a linear function of the speed = |v|, with a and b being constant.

Equation 7 of linear speed modulation is a reasonable approximation for place cells, at least when averaged over a population (McNaughton et al. 1983) (see also Fig. 8). A similar modulation by speed occurs for cells in the motor cortex of monkeys. Schwartz (1993, 1994) has shown that adding up the population vector head to tail approximately reproduces the hand trajectory. This implies that the length of the population vector is roughly proportional to reaching speed. Thus to a first approximation, the reaching speed modulates the average firing rate linearly also in the motor cortex.

View larger version (14K): [in this window] [in a new window]   FIG. 8. Modulation of place cell firing rate by running speed was nearly linear. A: histogram of speed distribution. Speed was computed after convolving position trajectory with a Gaussian kernel. B: average firing rate across all cells as a function of binned running speed. For animal 1, mean speed was 19.6 cm/s, mean firing rate was 0.92 Hz, and linear fit was 0.028 + 0.39 with correlation coefficient r = 0.991. For animal 2, corresponding parameters were 9.3 cm/s, 1.1 Hz, and 0.037 + 0.73, with r = 0.992. Only 1st half of data are shown (7 min for animal 1 and 10 min for animal 2). For 2nd half of data, parameters for animal 1 were 18.0 cm/s, 0.91 Hz, and 0.033 + 0.31 with r = 0.991; for animal 2, they were 7.0 cm/s, 1.38 Hz, and 0.023 + 1.25, with r = 0.956.

THETA RHYTHM. The timing of spike firing relative to the phase of the theta waves in the hippocampal local field potential encodes additional information about the location of the rat (Brown et al. 1996; Jensen and Lisman 1996; O'Keefe and Recce 1993; Samsonovich and McNaughton 1996; Skaggs et al. 1996; Tsodyks et al. 1996). Using phase information would be equivalent to splitting a single place field into several smaller regions while preserving the total number of spikes. How this might affect reconstruction accuracy is described quantitatively by Eq. 53. However, information from field potential recordings was not used here because with 20-30 simultaneously recorded cells, the optimal time window for reconstructing position was ~1 s (Fig. 6), which spanned around 8 theta cycles and averaged out the phase information from a single cycle.

View larger version (18K): [in this window] [in a new window]   FIG. 6. Mean error of Bayesian 2-step reconstruction method depends on time window and shift of its center. A: length of time window was altered systematically, whereas all other parameters were identical to those in Fig. 5 with all cells. B: center of time window (length fixed at 1 s) was shifted relative to current instant of time. Only data from first half session were used for sampling. Reconstruction was performed on data from both first half session (light gray curves and symbols) and second half session (black curves and symbols, with two outlier data points removed). For animal 1, each half session lasted 7 min and for animal 2 each lasted 10 min. For both animals, reconstruction for second half session was most accurate when center of time window was slightly shifted towards past.

	LINEAR APPROACH: BASIS FUNCTIONS

Direct basis and reciprocal basis

A basis function method uses a linear combination of fixed templates, or basis functions, with the coefficients in the linear combination proportional to the activity of the neurons. The template functions can be chosen arbitrarily. Define n_i as the number of spikes fired by cell i within the time window, and _i(x) as an arbitrary basis function or template function associated with this cell. The basic computation is the linear sum

(8)

which is a distribution function over two-dimensional (2-D) physical space indexed by x. The peak position of this distribution can be taken as the reconstructed position of the animal; that is

(9)

where argmax means the value of x that maximizes the function.

View larger version (41K): [in this window] [in a new window]   FIG. 2. Snapshots of reconstructed distribution density for position of a rat (animal 1) as compared with its true position. True position occupied a single pixel on 64 × 64 grid, corresponding to 111 × 111 cm in real space. Probabilistic method (Bayesian 1-step) often yielded a sharp distribution with a single peak, whereas direct basis method typically led to a broader distribution with multiple peaks. Total numbers of spikes collected from all 25 cells for three snapshots were 5, 22, and 31, respectively. Time window for each snapshot was 0.5 s.

The basis function framework encompasses existing methods such as the population vector method (Georgopoulos et al. 1986), the template matching method (Lehky and Sejnowski 1990; Wilson and McNaughton 1993), and the optimal linear estimator (Salinas and Abbott 1994) as special cases. A more detailed discussion of these methods follows in a later section.

The major question in this approach is how to choose template functions with minimal reconstruction error. It is useful to make a distinction between a direct basis and a reciprocal basis. A direct basis is identified with the firing rate distribution, which can be measured experimentally. A reciprocal basis, which arises naturally in minimizing the mean square error, is related to the Moore-Penrose pseudoinverse. Here reciprocal means that applying the pseudoinverse twice reverts back to the original basis.

To formalize the method, let f_i(x) be the average firing rate of the cell i when the animal is at position x. Then, by definition, the expected number of spikes n_i this cell fires within a time window of unit length should be

(10)

where the distribution function (x) describes the position of the animal and is 1 at the animal's current position and is zero elsewhere. Equation 10 remains valid even when (x) is an arbitrary probability distribution density.

The goal is to recover the position distribution (x) by combining some fixed template function _i(x) associated with each cell i. For an ideal reconstruction, we should have

(11)

for an arbitrary position distribution (x). In practice, we seek an approximate reconstruction by minimizing the error.

To find the optimal template functions g_i(x), we use vector-matrix notation to rewrite Eq. 10 as

(12)

where n = (n₁, n₂, . . . , n_N)^T is a column vector, F = [f₁, f₂, . . . , f_N] is a matrix with each column vector f_i obtained by concatenating all the pixels in the firing rate map f_i(x). The exact method of the concatenation does not affect the final result, as long as it is used consistently. I_k is the kth column of the identity matrix I. The sole nonzero element of I_k represents the pixel of the animal's position after the concatenation. The requirement Eq. 11 for perfect reconstruction now becomes

(13)

where, in the last step, Eq. 12 has been used, and matrixG = [g₁, g₂, . . . , g_N] is constructed in the same way as F, and each column vector g_i corresponds to the basis function g_i(x) to be determined. If Eq. 13 holds true for all k, that is, for all positions of the animal, then we must have I = GF^TI. Thus in general we should try to find G that minimizes

(14)

where X² means the sum of the squares of all the elements of matrix X. The solution to this least-mean-square problem is well known

(15)

where  denotes Moore-Penrose pseudoinverse and T denotes transpose. F and G are reciprocal in the sense that we also have

(16)

In other words, applying the pseudoinverse and transpose twice returns to the original basis. An example of reciprocal basis functions for place fields is shown in Fig. 1.

The reciprocal nature of the pseudoinverse can be more clearly seen using singular value decomposition (Golub and Van Loan 1996). Starting from the expansion F = UDV^T, where U and V are orthogonal matrices and D is a rectangular diagonal matrix, we immediately obtain G = UDV^T, where D is the diagonal matrix obtained from D^T by replacing each nonzero element _i with its inverse 1/_i, hence the reciprocal property. The familiar method of computing pseudoinverse from the normal equation

(17)

by inverting the correlation matrix F^TF yields identical results, provided that the matrix is not singular or near singular.

It follows from Eq. 17 that the firing rate maps themselves become identical to the optimal template functions if they do not overlap with each other. That is, F = G or f_i(x) = g_i(x) if

(18)

or equivalently, the tuning functions f₁(x), . . . , f_N(x) are uncorrelated

(19)

where C is a proportional constant and _ij is Kronecker delta that equals 1 when i = j and 0 otherwise. The conditions for orthogonality Eq. 18 or Eq. 19 are approximately satisfied for the recordings in Fig. 1, which explains why the reciprocal fields look similar to the place fields (direct basis) and why the corresponding performances are also similar (Figs. 3 and 5).

View larger version (62K): [in this window] [in a new window]   FIG. 3. True X and Y positions of animal 1 running on figure-8 maze as compared with positions reconstructed by different methods with 25 place cells. Same 60-s segment is shown in all plots. Time window for reconstruction was 0.5 s, which was moved forward with a time step of 0.25 s. For a fair comparison of different methods, if none of 25 cells fired within time window, reconstructed position at preceding time step was used. Probabilistic or Bayesian methods were especially accurate and erratic jumps in reconstructed trajectory were reduced by a continuity constraint by using information from two consecutive time steps.

View larger version (28K): [in this window] [in a new window]   FIG. 5. Mean errors of different reconstruction methods applied to data from two different rats. Errors all decreased as more cells were used, and lowest errors were comparable with intrinsic error of tracking system, which was about 5 cm (horizontal lines). Shaded region: reconstruction errors prohibited by lower bound derived by Fisher information and Cramér-Rao inequality (Eq. 53). Each data point and error bar represent average and standard deviation of mean errors in 40 repetitive trials in which a subset of cells were drawn randomly from whole sample. Animal 1 ran on figure-8 maze; 7 min of data were used for sampling and subsequent 7 min of data were used for reconstruction. Animal 2 ran on rectangular maze (upper loop of figure-8 maze) and both sampling and reconstructing time intervals were 10 min. Time window was 1 s.

In the above analysis, only the case where the animal is equally likely to visit any position x has been considered; that is, the probability distribution P(x) is assumed to be uniform. In this case, the direct basis is the firing rate map _i(x) = f_i(x), and the optimal template function that minimizes the mean square error is the reciprocal basis _i(x) = g_i(x) given by Eq. 15. For nonuniform P(x), instead of using the raw firing rate maps f_i(x) as the basis functions, choose

(20)

This weights the basis functions by the prior probability for the animal to visit each position. As expected, simulations confirm that this modification reduces reconstruction error. Similarly, we can use

(21)

as the reciprocal basis functions. Other modifications of the reciprocal basis are considered in APPENDIX A and compared in Table 1. In our comparison of different methods in Figs. 3 and 5, the simpler Eqs. 20 and 21 were used.

In summary, given the numbers of spikes (n₁, n₂, . . . , n_N) fired by all the cells during a given time interval, combine the chosen basis functions _i(x) linearly as n_ii(x), and then take the peak position x ^as the reconstructed position of the rat as described by Eq. 9. The whole trajectory is constructed by sliding the time window forward (Fig. 3).

Simple example of reciprocal basis

The reciprocal basis method is illustrated by a hypothetical example of coding in the motor system (Fig. 4). Imagine two motoneurons each of which, when activated alone, will drive movement in the direction f₁ or f₂, respectively. We call f_i the driving direction of neuron i. Assume that the two neurons are activated randomly and independently and that the actual movement is the vector sum

(22)

where n₁ and n₂ are the activity of the two neurons relative to their resting levels.

View larger version (16K): [in this window] [in a new window]   FIG. 4. A hypothetical two-neuron motor system showing distinction between driving directions and preferred directions, which form a pair of reciprocal vector bases. Apparent preferred directions differ from true driving directions unless system is orthogonal. Here reciprocity means that if vectors g1 and g2 are taken as driving directions, then vectors f1 and f2 will become apparent preferred directions.

Suppose that both the activity of one neuron and the overall movement are being measured. How is the activity of this neuron correlated with the actual movement direction? Each neuron i will appear to have a cosine directional tuning with a fixed preferred direction g_i, which, however, is different from the true driving direction f_i. To see why this is the case, write the driving directions f₁ and f₂ as column vectors and define F = [f₁, f₂] as a 2 × 2 matrix. Equation 22 is equivalent to

(23)

(24)

(25)

is the transposed inverse. This is a special case of the reciprocal basis in Eq. 15 because the pseudoinverse is identical to the ordinary inverse for an invertible square matrix. Expressing G = [g₁, g₂] in terms of two column vectors g₁ and g₂, the solution to Eq. 24 becomes

(26)

Therefore, the activity n_i of each neuron is the inner product between the actual movement direction f and a fixed direction g_i, which is the preferred direction. The driving directions and the preferred directions are reciprocal in the same sense as in the preceding section. They are identical only when the two driving directions are orthogonal. In general, Eq. 25 means F^TG = I (unit matrix), or equivalently, f₁·g₁ = f₂·g₂ = 1 and f₁·g₂ = f₂·g₁ = 0. That is, in general, the preferred direction g₁ of neuron 1 is always perpendicular to the driving direction f₂ of neuron 2, and vice versa. The distinction between a pair of reciprocal vector bases is equivalent to that between contravariant and covariant tensors (Pellionisz 1984; Pellionisz and Llinas 1985).

The biological meaning of a reciprocal basis is not always as clear as in the simple example given above. The reciprocal field for each place cell depends not just on the place field of that cell but also on all other cells in the population. If cells work together in groups or modules that are relatively independent of each other, then the reciprocal pairs should only include those cells in the same group. A different grouping can affect the shapes of the reciprocal fields. The reciprocal fields as shown in Fig. 1 are meaningful only with respect to those 25 cells in question and are therefore only a theoretical construct. Nonetheless, the concept of reciprocity they illustrate might have biological implications.

Special cases of basis functions

TEMPLATE MATCHING. The template matching method for reconstruction has been used for modeling stereo hyperacuity (Lehky and Sejnowski 1990) and place cells (Fenton and Muller 1997; Wilson and McNaughton 1993). This method is equivalent to the direct basis method. At each fixed position x, the average firing rates of all the cells can be described by a vector f = (f₁(x), f₂(x), . . . , f_N(x)). A vector template is stored for each position x and is later used to match the actual number of spikes fired by these cells within a time window: n = (n₁, n₂, . . . , n_N). During reconstruction, the matching method searches for the position in which the dot product n·f is maximized. This is equivalent to maximizing n·f = n_if_i(x) using f_i(x) as basis functions. In other words, the reconstructed position is

(27)

Uniform scaling by the numbers of spikes does not affect the peak position or the reconstruction result. Position-dependent scaling such s(x) = 1/f_i(x)² can be accommodated by redefining each template function f_i(x) as f_i(x)s(x).

VECTOR METHODS. Vector reconstruction is widely applied to data from neural populations. Vector methods such as the population vector and the optimal linear estimator compute

(28)

where x_i is a fixed vector for each neuron i, and n_i is its activity level, such as the total number of spikes during a time interval or the firing rate.

This method is different from the basis function method because at each time step, the basis function method yields a scalar distribution function of the animal's probable position over the entire maze, as shown in Fig. 2. By contrast, in the same situation, a vector method would yield a single vector that specifies a single point on the maze. On the other hand, after being discretized and concatenated, the relation between the direct and reciprocal basis functions is similar to that between the population vector and the optimal linear estimator (Georgopoulos et al. 1986; Salinas and Abbott 1994; Sanger 1994).

The basis function framework includes the population vector method as a special case in the sense that the performance of population vector decoding is identical to combining cosine functions as the template functions. The direction of the population vector _vector in Eq. 28 is identical to the peak position of the following scalar function

(29)

where vector x is a free variable in the same space as _vector but has unit length, and the template function

(30)

is a cosine function of the angle between the two vectors. In other words, the direction of the reconstructed vector defined by Eq. 28 also can be computed by

(31)

where |_vector| is the length of the vector _vector. To see this, take the dot product of the both sides of Eq. 28 with the unit vector x and then use the definitions in Eqs. 29 and 30 to obtain (x) = x·_vector, which is maximized only when the unit vector x is aligned with _vector. Therefore basis functions with cosine tuning can precisely implement a vector method (see also Salinas and Abbott 1994).

POPULATION VECTOR WITH SCALING. For place cells, the population vector needs to be scaled by the total activity so that the reconstructed position

(32)

where x_i is the center of the place field of cell i. This method has been considered by several authors (Abbott and Blum 1996; Blum and Abbott 1996; Burgess et al. 1994; Muller et al. 1987). It will be shown later that the scaling in Eq. 32 can be justified by the Bayesian method.

The two forms of population vectors in Eqs. 28 and 32 have different physical meanings. A vector representing direction will be called a directional vector, and a vector representing spatial position will be called a positional vector. If we scale a directional vector by an arbitrary scalar, the resultant vector still points in the same direction, but scaling a positional vector will lead to a different position. Vector reconstruction methods are most suitable for reconstructing a directional vector. When reconstructing position in 2-D space, a vector method can produce implausible results, even after scaling by total activity (Fig. 3). We return to this point later.

To further illustrate the difference between the basis function method and a population vector, consider the special case in which every place field has the same Gaussian tuning function f_i(x) = G(x  x_i), where x_i is the center of the place field of cell i. To find the peak position of n_if_i(x), set its derivative with respect to x to zero. This leads to

(33)

which can be compared with Eq. 32. Because eachG(_basis  x_i) is small unless the estimated position _basis is sufficiently close to the place field center x_i, the overall effect may be interpreted as a population vector with an additional weighting scheme in favor of those place fields centered close to the estimated position. Consequently, an occasionally active cell with a place field far away from the estimated position does not affect the final result of the basis function method as much as in the case of a population vector.

	PROBABILISTIC APPROACH: BAYESIAN METHODS

Background

The Bayesian approach is a natural one for reconstruction. It is optimal within a probabilistic framework and, as shown here, yields the best reconstruction results for our dataset. Once the position dependence of the firing rates of a group of place cells is known, the Bayes formula directly addresses the inverse problem: given the firing rates of these cells, how can we infer the most probable position? Bayesian reconstruction has been applied by Földiák (1993) to visual orientation tuning, by Sanger (1996) to motor directional tuning, and by Ginzburg (Gerrard et al. 1995) and Brown (Brown et al. 1996) to place cells.

Two optional supporting models can be used with Bayesian reconstruction. The first is the tuning model, which is especially useful when data sampling is sparse. Brown, Frank, and Wilson (1996) used Gaussian tuning model in their reconstruction, and we used a Gaussian tuning model for estimating the information-theoretical limit on reconstruction accuracy. In our reconstruction, the tuning functions were determined empirically without assuming an explicit analytical model. A second option for Bayesian reconstruction is the spike generation model. Because of the low firing rates of place cells, we used the Poisson firing model in our reconstruction. Bayesian reconstruction using Poisson spike statistics was first used by Sanger (1996). Alternatively, without using an explicit spike generation model, the stimulus-response relation can be obtained empirically as a look-up table (Földiák 1993).

Basic method

Assume that using the methods in the preceding section we have measured the firing rate maps f₁(x), f₂(x), . . . , f_N(x) of a population of N place cells as well as the animal's position distribution P(x) during a period of time. At an arbitrary later time t, given the numbers of spikes fired by these place cells within the time interval from t  /2to t + /2, where  is the length of time window, the goal is to compute the probability distribution of the animal's position at time t. Notice that what is to be computed here is a distribution of positions instead of a single position; we always can take the most probable position, which corresponds to the peak of the probability distribution, as the most likely reconstructed position of the animal (Fig. 2).

Let the vector x = (x, y) be the position of the animal, and the vector n = (n₁, n₂, . . . , n_N) be the numbers of spikes fired by our recorded cells within the time window, where n_i is the number of spikes of cell i. The reconstruction is based on the standard formula of conditional probability

(34)

The goal is to compute P(x|n), the probability for the animal to be at position x given the numbers of spikes n. P(x) is the probability for the animal to be at position x, which can be measured experimentally. The probability P(n) for the numbers of spikes n to occur can be determined by normalizing the conditional probability P(x|n) over x and therefore does not have to be estimated directly.

The key step is to evaluate the conditional probability P(n|x), which is the probability for the numbers of spikes n to occur given that we know the animal is at location x. It is intuitively clear that this probability is determined by the firing rate maps of the place cells. More precisely, if we assume that the spikes have Poisson distributions and that different place cells are statistically independent of one another, then we can obtain the explicit expression

(35)

where f_i(x) is the average firing rate of cell i at position x, and  is the length of the time window.

The final formula obtained from Eq. 34 reads

(36)

The Bayesian reconstruction method uses Eq. 36 to compute the probability P(x|n) for the animal to be at each position x, given the numbers of spikes n of all the cells within the time window. In this probability distribution, the peak position or the most probable position is taken as the reconstructed position of the animal. In other words

(37)

By sliding the time window forward, the entire trajectory of the animal's movement can be reconstructed from the time varying activity in the neural population.

Testing the assumptions

The Bayesian reconstruction as applied in Eq. 36 relies on two assumptions, namely, Poisson spike distribution and independent firing of different cells. These assumptions are reasonable and convenient approximations, although neither is precisely true.

POISSON SPIKE STATISTICS. The Poisson firing model commonly is used for neuronal spike statistics (Tuckwell 1988). It is equivalent to assuming that the firing is as random as possible. If points are distributed randomly over a long time interval, with equal chance for each point to fall at any location independently of each other, then the number of points contained in a small time window has a Poisson distribution. A direct test for a Poisson distribution is not easy for place cells because their mean firing rates are low and depend on the location, which an animal may visit only transiently. Nonetheless, the interspike intervals can be compared directly with the exponential distribution expected for a Poisson process. Consider cell i with firing rate map f_i(x). According to the inhomogeneous Poisson model with mean firing rate depending on the position x, the overall probability density for interspike interval s should be

(38)

where P(x) is the probability for the animal to visit position x. The above formula implies that dp_i(s)/ds < 0 always holds true; that is, the interspike interval histogram must be monotonically decreasing.

Applying Eq. 38 to real place cell data, we found that it can approximately account for the distribution of intervals >10 ms. But there are several problems. First, the actual occurrence of briefer intervals was much more frequent than expected. This means that the Poisson model fails to capture the tendency for hippocampal pyramidal cells to fire in bursts, called complex spikes (Ranck 1973). Second, instead of decreasing monotonically, the actual interspike interval histograms tended to have a small bump within the theta frequency range (4-12 Hz), reflecting the periodic drive from the septal pacemakers (Stewart and Fox 1990). Third, the Poisson spiking model does not have a refractory period to rule out extremely short intervals. A related finding by Fenton and Muller (1997) is that place cell firing for similar movement trajectories in an open environment was more variable that expected from a Poisson distribution. Despite these limitations, the inhomogeneous Poisson model proved to be adequate for the purpose of reconstruction.

INDEPENDENCE OF DIFFERENT CELLS. The activities of two cells are statistically independent only if the joint probability of firing is a product of the individual probabilities. For two place cells, 1 and 2, their independence means that at each position x

(39)

where n₁ and n₂ are the number of spikes collected from the two cells within an arbitrary time window. For the data in Fig. 1, the majority of the pairs of place cells were approximately independent simply because most place fields did not overlap. In other words, most cell pairs seldom fired together. Thus Eq. 39 may hold for the trivial reason that both sides are zero. A related consequence is that multiplying place fields in Bayesian reconstruction typically leads an extremely sharp distribution, with zero value for most pixels (Fig. 2 and 13). The independence is no longer exactly true for cells with highly overlapping place fields (cf. Eichenbaum et al. 1989). In fact, the correlations in the firing of these cells could be strengthened rapidly by Hebbian-type learning process (Blum and Abbott 1996; Mehta and McNaughton 1997; Skaggs and McNaughton 1996; Wilson and McNaughton 1994). In such a situation, the independence assumption is only a simplifying approximation for theoretical convenience. There is additional information in correlated firing that could be exploited to improve the accuracy of reconstruction.

Continuity constraint: using two time steps

The trajectories reconstructed by all methods suffer from occasional erratic jumps, which often are caused by low instantaneous firing rates of the recorded place cells, especially when the animal stops running (Fig. 3). Clearly, if all recorded cells stop firing, there is not enough information for accurate reconstruction.

Introducing a continuity constraint can improve reconstruction accuracy by reducing erratic jumps in the reconstructed trajectory. The possibility that some of the jumps may reflect real biological processes will be discussed later. To implement the continuity constraint, the reconstructed position from the preceding time step is used as a guide for estimating the current position. Speed information can be used to limit the change in position allowed within a single time step.

The continuity constraint can be incorporated easily into both the Bayesian and the basis function frameworks. Here we demonstrate the constraint only for the Bayesian method. The one-step Bayesian method presented previously computes the conditional probability P(x_t|n_t) for the animal's position x_t at time step t, given the numbers of spikes n_t within a time window centered also at time t. Here subscripts denote the time step. When the reconstructed position x_t1 at the preceding time step is used as an additional piece of information, we can compute the conditional probability P(x_t|n_t, x_t1) from the general Bayes formula P(x_t1|x_t, n_t)P(x_t, n_t) = P(x_t|n_t, x_t1)P(n_t, x_t1) and the assumption that

(40)

In the sense of probability, this assumption means that the activity n_t at the current time step cannot directly affect the position x_t1 at the preceding time step. With the help of simple equalities like P(x_t, n_t) = P(x_t|n_t) P(n_t) and P(n_t, x_t1) = P(x_t1|n_t)P(n_t), we immediately obtain the final equation

(41)

where k = 1/P(x_t1|n_t) is a scaling factor that does not depend on x_t and can be determined by normalizing P(x_t|n_t, x_t1) over x_t.

Now the reconstructed position at time t is

(42)

Compared with the single step Eq. 36, the only new factor in Eq. 41 is the conditional probability P(x_t1|x_t). Intuitively, the continuity constraint should require that given the current location x_t, the preceding location x_t1 cannot be too far away. For a precise formulation, we assume a 2-D Gaussian distribution

(43)

where the standard deviation _t is related to the speed _t at time t by the formula

(44)

In general, if the standard deviation _t of the Gaussian prior is too large, the continuity constraint has little effect. On the other hand, if _t is too small, the constraint may become too restrictive and the reconstructed position might get stuck in the same position if the real position has moved away by a distance much larger than _t. In our analysis, we simply used d = 1 for linear movements and confined _t to between 20 and 60 cm (good empirical values). The scaling in Eq. 44 was chosen so that _t would just reach the allowed maximum of 60 cm at top running speed. For generality, we did not include any bias for the direction of movement. As shown in Figs. 3 and 5, the continuity constraint is effective in reducing reconstruction errors caused by erratic jumps.

Incorporating firing rate modulation

Under the assumption of multiplicative modulation of firing rate such as in Eq. 6, the Bayesian Eq. 36 should be modified by replacing the fixed tuning function f_i(x) of each cell i by

(45)

where, for simplicity, the modulation factor m(t) is supposed to be the same for all cells. Here m(t) is written as a function of the time step. This general formulation includes more specific variables such as modulation by speed as a special case. Note that the modulation factor can be estimated directly from spike data using the approximation.

(46)

where F(t) is the instantaneous firing rate averaged over all the cells and F_mean is the average of F(t) over time. The modified Bayesian formula reads

(47)

Relationship to basis functions

GENERAL EQUIVALENCE. In general, the performance of all basis function methods can be mimicked by the Bayesian method if appropriate template functions are used. To implement the sum n_i(x) with the basis functions _i(x), choose the tuning function

(48)

and the prior probability for spatial occupancy

(49)

(50)

Conversely, the Bayesian method also can be implemented in a basis function framework, provided that a constant bias is allowed. This relationship will be considered in detail in our later discussion on network implementations.

RELATIONSHIP WITH POPULATION VECTOR WITH SCALING. Under special assumptions, the Bayesian method can lead to the population vector method with scaling by total activity. This scaling is needed for reconstruction of a positional vector. Suppose the tuning function for each cell i is a Gaussian f_i(x) = A_i exp(x  x_i²/2²_i), where x_i is the center of the place field of cell i and _i characterizes its width. Suppose the centers x_i and the peak firing rates A_i are uniformly distributed in space, then the sum f_i(x) is approximately constant in space. In this case, the optimal position in Bayesian reconstruction is given by

(51)

where the second equality can be obtained by setting the derivative of the product with respect to x to zero. This is a generalized population vector, which reduces to the form in Eq. 32 when all place fields have identical width; that is, when all _i = . This justifies scaling the positional population vector by the total activity n_i in Eq. 32.

	RESULTS OF PLACE CELL RECONSTRUCTION AND DISCUSSION

Comparison of different methods

RECONSTRUCTION ACCURACY. The performance of different methods can be compared qualitatively in Fig. 3 and quantitatively in Table 1 and Fig. 5. The effects of various parameters are discussed later.

The reconstruction error of the direct basis method for a given number of cells is similar to that reported by Wilson and McNaughton (1993). This result makes sense because of the essential equivalence between the matching method and the direct basis method as discussed before.

The linear method with a reciprocal basis only improved slightly the performance of the direct basis because the place cells were mostly orthogonal and the two sets of basis functions were quite similar to each other (Fig. 1). Further modifications of the reciprocal basis sometimes brought additional improvement (Table 1).

Probabilistic methods based on the Bayesian approach were more accurate. The best method overall was the Bayesian method using two time steps to enforce a continuity constraint on the reconstructed position. Its errors were in the range of the intrinsic error of the system for tracking the animal's position, which was estimated to be ~5 cm (Skaggs et al. 1996; Wilson and McNaughton 1993).

Because reconstruction errors tend to be much larger during immobility (see further text), excluding slow running periods from the analysis can improve reconstruction accuracy (Gerrard et al. 1995). However, data exclusion requires additional criteria. For all methods in this paper, we reconstructed the whole trajectory without introducing additional parameters.

PROBLEM WITH VECTOR METHODS. The population vector method was the least accurate in reconstructing position in space for our dataset. The major problem is that it sometimes yielded implausible results. To see why, imagine that only two cells were active within the time window; then the weighted average of their positions would lie somewhere on the line connecting the two centers of the place fields. The estimated position may be in the center of the maze, which was never visited by the rat. This is a general problem for all vector methods including the reciprocal vector or optimal linear estimator. As long as the shape of the environment is not convex, a vector method may suffer from this problem. The existence of systematic bias in vector methods has been pointed out by Muller, Kubie, and Ranck (1987). In contrast, both the basis function method and the probabilistic method take into account the entire 2-D distribution functions, and the final peak always occurs in a plausible position. See also the discussion after Eq. 33.

REMARK ON SILENT TIME STEPS. For a fair comparison of different methods, in all of our simulations we have used the heuristic that if all the recorded cells were silent within a time window, the current reconstructed position was set to be the reconstructed position at the preceding time step. We needed this rule because when all recorded cells are silent within the time window, all basis function methods fail, whereas the Bayesian methods are robust. With a time window of 0.5 s, ~15% of the time steps contained no spike in animal 1, and the percentage was only ~0.2% in animal 2. For a time window of 1 s, the percentages dropped to 4.3 and 0% for animals 1 and 2, respectively. There were fewer silent time steps in animal 2 because its maze was smaller, but there were more cells so that the overall average firing rates were higher than that of animal 1.

Effects of parameters

NUMBER OF CELLS. When more cells were used in reconstruction, the errors tended to decrease (Fig. 5). Different methods maintained their relative accuracy for different numbers of cells. In theory, the error of various reconstruction methods should be inversely proportional to the square root of the number of cells (Eq. 56). However, the inverse square root law is valid only when a sufficiently large number of cells is used, as shown by simulation (Fig. 12). Here the number of cells was too small to apply the law.¹

View larger version (36K): [in this window] [in a new window]   FIG. 12. Comparison of errors of population vector and Bayesian methods for 2-D reaching problem assuming cosine tuning, Poisson spike statistics, and a uniform probability distribution of preferred directions around a circle. Simulation results (symbols) agreed with theoretical results (curves and straight lines in log-log plot, given by Eqs. 66 and 68) when number of cells N was large so that errors of both methods dropped inversely proportional to . In particular, Cramér-Rao lower bound (lower theoretical curve or line) was reached asymptotically by Bayesian method when there were more than 20 cells. Parameters: fmax = 10 Hz, fmin = 1 Hz and each data point is average of 104 trials with randomly chosen preferred directions.

TIME WINDOW. A time window with a duration of ~1 s was a good choice for the Bayesian two-step method for these data (Fig. 6A). The results of other methods were similar (not shown), although the errors for large time windows sometimes leveled off or only increased slightly. In theory, the error is expected to be inversely proportional to the square root of the time window (Eq. 56). This can account for the general tendency for the reconstruction to improve with a longer time window. However, this theoretical consideration no longer holds when the time window exceeds the optimal size.

TIME SHIFT. In the Bayesian two-step method, the center of the time window was shifted relative to the current time to see what shift produced the best accuracy (cf. Blair and Sharp 1995; Blair et al. 1997; Muller and Kubie 1989; Taube and Muller 1995). Alignment of the time window to within 100 ms of the current time gave the best results. Figure 6B compares the reconstruction errors from the first half of the session to that from the second half. As expected, the errors for reconstructing the second half were larger than that of self-reconstruction because all parameters for the reconstruction algorithm were determined using data from only the first half session.

Upon closer examination of Fig. 6B, it would appear that the minimum of the curve for reconstructing the second half session (dark symbols) drifted toward the past with respect to the minimum of the self-reconstruction curve (light symbols). To quantify the drift, each curve was fitted with a fifth-order polynomial, and the minimum of the smooth curve was found. The drift was computed as the difference between the two minima. For the data shown in Fig. 6B, the value of the drift was 66 ms (34  32 ms) for animal 1 and 68 ms [93  (25) ms] for animal 2. Conversely, when the second half of data was used to reconstruct the trajectory of the first half, the minimum appeared to drift towards the future with respect to the minimum of self-reconstruction (not shown). Using the same smoothing method, the value of the drift was 118 ms [107  (11) ms] for animal 1 and 93 ms [17  (76) ms] for animal 2. 

View larger version (26K): [in this window] [in a new window]   FIG. 7. Mean reconstruction errors of Bayesian 2-step method in consecutive and nonoverlapping time blocks, each lasting 3 min. Only a single block (shaded) was used as data block for sampling and reconstruction was performed on all blocks. Graphs for two animals are aligned along data blocks. Average reconstruction errors tended to change gradually even though there was no overlap between contiguous time blocks.

The accuracy for the estimation of the drift was limited by noise and the flatness of the curves. Nevertheless, the direction and the magnitude of the drift were consistent with the phenomenon of the back shift of place field centers during repetitive running as reported by Mehta and McNaughton (1997), who directly computed the center of mass of place fields and found an average back shift of a few centimeters during a running period of the order of several minutes. On the other hand, we compared the spatial distribution of the reconstruction errors as a vector field for different time blocks (see further), but could not see an obvious systematic shift of the error direction in correlation with the running direction.

TIME BLOCKS: SLOW TRENDS. As shown in Fig. 7, when a single 3-min time block of data was used to reconstruct the movement trajectory of other time blocks, the average reconstruction error tended to drift systematically. That is, reconstruction tended to get worse when going further into either the past or the future. Because there was no overlap between the time blocks, the tendency for gradual drift of error must be caused by gradual changes of either behavioral or neural origin.

During repetitive running, the mean running speed of both animals slowly decreased. For animal 2 (on the rectangular maze), the mean speed dropped from 9.0 cm/s in the first 3-min time block to 5.9 cm/s in the last time block, with an overall average of 8.2 cm/s across all seven time blocks. Similarly, the speed for animal 1 (on the figure-8 maze)dropped from 21.8 to 16.6 cm/s, for an average of 18.9 cm/s.On the other hand, the mean firing rate increased almost monotonically for animal 2 (from 0.86 to 1.34 Hz, for an average of 1.23 Hz) but showed less clear trend for animal 1 (from 0.94 to 0.87 Hz, for an average of 0.92 Hz). Therefore, slow changes both in behaviors and in neural properties could contribute to the gradual error drift in Fig. 7.

FIRING RATE MODULATION. Including firing rate modulation in reconstruction had only a small effect on the accuracy, even though the average firing rate was modulated clearly by running speed (Fig. 8), as reported by McNaughton, Barnes, and O'Keefe (1983).

Including firing rate modulation has no effect on the performance of basis function methods and the population vector, provided that an identical modulation factor is used for all cells as in Eq. 45. But under the same conditions, the performance of the modulated Bayesian method (Eqs. 46 and 47) may slightly improve (Table 1).

OTHER PARAMETERS. The reconstruction results were not sensitive to parameters like the exact width of the spatial grid and the exact radius of the blurring kernel as long as they were kept below the intrinsic error of the tracking system. In all reconstruction data presented in this paper, we used64 × 64 grid, corresponding to 111 × 111 cm in real space. Adding a small constant background value to the firing rate maps sometimes slightly improved the Bayesian reconstruction.

Finally, in the Bayesian two-step method, the standard deviation _t of allowable jumps in Eq. 44 varied between 20 and 60 cm, which gave good performance for the data of animal 1. Smaller values slightly improved the reconstruction for animal 2, which did not have many erratic jumps in the first place. Changing these values affects the reconstruction error only slightly, as long as they are not too small or too large.

Discontinuity in reconstructed trajectory

As shown in Figs. 3 and 9, reconstruction errors were relatively small most of the time but were punctuated by large jumps. Most of these jumps occurred when the animal stopped running (Fig. 9). This also can be seen in Fig. 3, where the jumps tend to occur between the corners, the places where the animal stopped for food reward. The results from animal 2 were qualitatively similar, but the jumps in reconstruction errors were less frequent because more cells were used in reconstruction and their average firing rates were higher.

The occurrence of these jumps is correlated with the several confounding factors, all of which are apparent in Fig. 9: 1) lower running speed or immobility, which itself is known to be correlated with; 2) lower firing rates of place cells (McNaughton et al. 1983) (see also Figure 8); 3) occurrence of large amplitude irregular activity and more frequent sharp waves in local field potentials (Buzsáki 1986; Vanderwolf et al. 1975) (in Fig. 9, the sharp waves can be seen as the thin vertical dark bands ~150-200 Hz in the Fourier power spectra of channel 1 of the local field potentials); and 4) the disappearance of the theta rhythm (Vanderwolf et al. 1975) (in Fig. 9, theta modulation can be seen as the fuzzy horizontal dark bands ~8 Hz in the Fourier power spectra of channel 2 of the local field potentials).

Because of the correlated factors listed above, reconstruction tended to be more accurate during the theta rhythm when the animal was running (Gerrard et al. 1995). This is consistent with the report that place cell firing represented an animal's position somewhat more faithfully during the theta rhythm (Kubie et al. 1984).

Another related factor is the sharp waves, which are believed to originate in the CA3 region perhaps because of the lateral connectivity in this region (Buzsáki 1986, 1989). Sharp waves are more frequent during immobility. It is possible that the activity during sharp waves might be involved with other functions such as memory consolidation and reactivation of learned correlations rather than merely reflect the current spatial position of the animal. Reconstruction based on sharp wave activity may yield an incorrect position. However, because the jumps did not always coincide with the occurrence of sharp waves, we can only conclude that sharp waves are one possible contributing factor for discontinuity in reconstructed path.

A further consideration concerns the hippocampal view cells found in monkey (Ono et al. 1993; Rolls et al. 1995). One cannot completely rule out the possibility that the activity of some cells in rat hippocampus might be similarly modulated by the animal's view rather than by the current position alone. If this is true, the reconstructed position may jump when a rat looks around or shifts attention.

In summary, most erratic error jumps occurred while the animal was not running. Many of these jumps can be accounted for by a momentary drop of instantaneous firing rates of recorded cells because fewer spikes imply less information for reconstruction. Some other jumps might reflect real biological activity in the hippocampus when the animal stopped running to eat, look around, and plan the next move.

	RECONSTRUCTION ACCURACY: THEORETICAL LIMITS AND BIOLOGICAL IMPLICATIONS

Fisher information

The accuracy of different reconstruction methods varied greatly. What is the maximum accuracy that can possibly be achieved? The problem is to estimate the maximum amount of information inherent in the neuronal spikes that can be extracted for reconstruction.

Fisher information is particularly suitable for addressing this question (Paradiso 1988; Seung and Sompolinsky 1993; Snippe 1996). It differs from Shannon information, which has been used more commonly in estimating information contained in spikes (for a recent review, see Rieke et al. 1997). Although these two information measures are related, their relationship is not straightforward (Cover and Thomas 1991). Fisher information depends on the shape and the slope of the tuning function (see APPENDIX B). In the place cell example, if the spatial bins are scrambled, the amount of Fisher information would be different. By contrast, Shannon information would be the same regardless of the ordering because it is not directly related to the slope of a distribution function. Shannon information offers a useful measure on the distribution of place fields (Skaggs et al. 1993) but no existing theory has linked it to reconstruction accuracy. See also Treves, Barnes, and Rolls (1996) for related discussion.

A key property of Fisher information is that its inverse is a lower bound, called the Cramér-Rao bound, on the variance of all unbiased estimators (Cramér 1946; Kay 1993; Rao 1965; Scharf 1991; Zacks 1971). This directly yields an estimate of the minimal reconstruction error that can possibly be achieved.

Minimal error under Gaussian tuning in 2-D space

To estimate Fisher information in place cell data, several properties of the place fields must be provided: the average shapes of place fields, the spiking statistics of the cells, the density of the place fields per unit area, and the statistics of the tuning parameters. Instead of using noisy real data for direct numerical evaluation, we use analytical models based on simplifying assumptions. Although these assumptions are not rigorously true, they are reasonably good approximations that allow useful analytical estimates to be derived. First, assume that the spatial tuning function is a 2-D Gaussian

(52)

Here we only state the final results of the detailed derivation provided in APPENDIX B. The simplest result is that reconstruction error cannot be smaller than

(53)

where ₂ = /2 is a constant correction factor, ² is the average of the square of the tuning width for different place fields, and N_spikes is total number of spikes collected from all the cells within the time window, namely,

(54)

Here N is the total number of cells recorded, f_mean is their mean firing rate, and  is the length of the time window. The mean firing rate is related to the parameters of the Gaussian place fields by the formula

(55)

where A is the area of interest, in which the N place fields are uniformly distributed, and f_max is the average peak firing rate of all place cells. Here the assumption of the independence of peak firing rate and the tuning width has been used.

Another equivalent error formula is

(56)

The equivalence of the two Eqs. 56 and 53 follows immediately from Eqs. 54 and 55.

Equations 53 and 56 are expected to be more reliable when a larger number of cells is used. This is because when the sample size is too small, a reconstruction method may suffer from systematic bias; that is, for a given true position, the average reconstructed position in repetitive trials does not approach the true position. In addition, in our derivation of the error formulas, the procedure of replacing a sum by an integral makes sense only when a large number of cells is used (APPENDIX B).

Unbiased reconstruction is assumed implicitly when the Cramér-Rao bound for variance is used to estimate the minimal reconstruction error. Otherwise, it would be trivial for a reconstruction method to achieve zero variance simply by giving a constant output regardless of the input. Because we used no more than 30 place cells in the reconstruction, systematic bias did exist for all reconstruction methods. This factor may reduce the accuracy of the error formulas (cf. Fig. 12).

In deriving these formulas, we have assumed implicitly that the duration of the time window is small compared with the running speed so that the rat does not move too far during that period of time. So if the time window is too large compared the running speed, the formulas are no longer reliable.

Minimal error under Gaussian tuning in arbitrary dimensional space

To generalize the results in the preceding section, we can derive error formulas for an arbitrary spatial dimension D by assuming Gaussian tuning function, Poisson spike statistics, uniform distribution of the center of the Gaussians in space, and statistical independence of the firing of different cells as well as their peak rate f_max and tuning width . As the counterpart of Eq. 53, the minimal reconstruction error based on Fisher information is

(57)

where   means average over all the cells. The counterpart of Eq. 56 is

(58)

where constant coefficient

(59)

depends only on the dimension D. Here, as before, _D is the correction factor for dimension D (see Table B1 and Eq. B26);  is the width and f_max is the peak firing rate for the Gaussian tuning function (52) in D-dimensional space; and N_spikes is the total number of spikes collected from N cells within the time window

(60)

The mean firing rate f_mean averaged over all cells is related to the tuning parameters by

(61)

The equivalence of Eqs. 57 and 58 follows from Eqs. 60 and 61. The new parameter  is the cell density. For example, suppose N cells are distributed randomly within a D-dimensional cube with edge length L, we have

(62)

In the 1-D case,  = N/L is the number of cells per unit length; in the 2-D case,  = N/L² is the number of cells per unit area; and in the 3-D case, = N/L³ is the number of cells per unit volume. Parameter  is a characteristic length so that ^D is the average length, area, or volume per cell when D = 1, 2, or 3, respectively (Fig. 10A). When D = 2; the new Eqs. 57 and 58 reduce to the old Eqs. 53 and 56 for place fields.

View larger version (18K): [in this window] [in a new window]   FIG. 10. A: distribution of centers of Gaussian tuning functions in 1-, 2-, and 3-D spaces. Only regular grids are shown here. In actual simulation, centers were scattered randomly with a uniform distribution. Average length, area, and volume per cell are given by , 2, and 3 for 1-, 2-, and 3-D cases, respectively. B: minimal achievable reconstruction error could increase, decrease or stay same as width of Gaussian tuning function increases, depending on dimension of space. Theoretical curves (Eq. 58) are compared with numerical simulation (symbols) of Bayesian reconstruction assuming model cells with Poisson spike statistics and identical Gaussian tuning function. Each data point is average of 104, 103, and 102 repetitive trials for 1-, 2- and 3-D cases, respectively. Without loss of generality, we set  = 1 in all simulation. Time window  = 1 s and peak firing rate fmax = 10 Hz.

View larger version (22K): [in this window] [in a new window]   FIG. 11. Mean angular errors for different test directions in a simple system with only 2 fixed orthogonal basis vectors, indicated by arrows in polar plots. A perfect cosine tuning curve and Poisson spike statistics were assumed. Each data point is an average over 5 × 104 repetitive trials. In this example, population vector was more accurate than Bayesian method when test direction was close to 225°, but its average error across all directions (26.47°) was still larger than that of Bayesian method (25.43°). Parameters: fmax = 10 Hz, fmin = 1 Hz, and time window  = 1 s.

The minimal reconstruction error can be reached by the Bayesian method (cf. APPENDIX C). In Fig. 10, the theoretical curves are compared with the simulation results of Bayesian reconstruction, assuming model cells with Poisson spike statistics and Gaussian tuning functions with randomly distributed centers but identical tuning width and peak firing rate. Here the error of the Bayesian method practically has reached the theoretical lower bound as long as the tuning width  was not too small compared with the characteristic length  = 1. In this example, it follows from Eq. 58 that

(63)

In other words, as the tuning width  increases, the error increases with in the 1-D case, keeps constant in the2-D case, and decreases with 1/ in the 3-D case.

Minimal error under cosine tuning and comparison with population vector

Cosine tuning is widely used to model motor cortical cells (Georgopoulos et al. 1988). We will show that under cosine tuning, the Bayesian method can achieve the best possible accuracy by reaching the theoretical lower bound, whereas the performance of the population vector is well above the bound.

However, the population vector can be more accurate than the Bayesian method in special situations. An example is given in Fig. 11, where two cells had orthogonal preferred directions, perfect cosine tuning, and Poisson spike statistics. The population vector outperformed the Bayesian method in a small region where the test direction was the furthest away from the two preferred directions.

Consider the 1-D situation that corresponds to a 2-D work space for reaching. We have N idealized cells with identical cosine tuning curves so that the mean firing rate of cell i is described by

(64)

where  is the angle for the test direction and _i is the preferred direction, which is drawn randomly from a uniform distribution around the circle, and

(65)

are constants. Suppose the spikes have Poisson statistics and different cells are independent. The average reconstruction error _Bayes for the Bayesian method Eq. 36 is expected to approach the minimal achievable error defined by Fisher information and Cramér-Rao bound, provided that a large number of cells is used (APPENDIX C). The minimal angular error (arc) derived by Fisher information is

(66)

where ₁ is the 1-D correction factor (Table B1), N is the total number of cells, and

(67)

In the same situation, the average angular error of the population vector method also can be estimated directly (APPENDIX D). The result is

(68)

where

(69)

with A and B given by Eq. 65. In the 2-D case, corresponding to 3-D work space for reaching, Eq. 68 needs to be modified by replacing ₁ by ₂ and Q₁ by

(70)

Because the error for the population vector in Eq. 68 and the error for the Bayesian method in Eq. 66 are both inversely proportional to the square root of the number of cells (see also footnote 1), the ratio

(71)

is independent of the number of cells. It can be verified analytically that J₁Q₁ < 1 regardless of the values of f_min, f_max, and the time window . This means that here the Bayesian method is always more accurate than the population vector regardless of the choice of parameters.

Place cell reconstruction error compared with theoretical limits

COMPARISON WITH LOWER BOUND. As shown in Table 1 and Fig. 5, our best reconstruction error by the Bayesian two-step method using all the cells was larger than the theoretical limit in Eq. 53 by a factor of about 2. To compute the theoretical lower bound, the width of the place fields was estimated using the first half of data by the method described in APPENDIX B. For animal 1, the average width of place fields was  = 10.1 cm and  = 11.2 cm. For animal 2,  = 8.6 cm and  = 9.6 cm. The mean firing rates of the 25 cells from animal 1 was f_mean = 0.92 Hz; and the average of the 30 cells from animal 2 wasf_mean = 1.09 Hz.

Because the Bayesian two-step method uses information from two time steps, is it fair to compare its performance with the Fisher information limit evaluated with the time window of a single time step? This still may be a reasonable comparison because the information from the previous time step is very crude: in actual reconstruction the width of the Gaussian prior distribution varied from 20 to 60 cm. It improved reconstruction accuracy mainly by prohibiting erratic jumps caused by low firing rates rather than by providing precise location information. For example, it can be seen in Fig. 3 that the Bayesian one-step method, while suffering from the erratic jumps, already provides accurate information about the fine details of the movement trajectory. In animal 2, the average firing rate for which was higher, the continuity constraint from the preceding time step actually did not contribute much to the accuracy(Fig. 5).

In idealized situations, the performance of the Bayesian method should reach the Cramér-Rao lower bound (Figs. 10 and 12 and APPENDIX C). For real place cell data, the basic assumptions of Poisson spike statistics and independence for the Bayesian formulas are not exactly true, which implies suboptimal performance. Reconstruction bias caused by the relatively small number of cells is another source of error (cf. Fig. 12). Finally, the estimate of Cramér-Rao bound itself relies on some simplifying assumptions such as Gaussian tuning, which may lead to additional error.

CORRECTION FACTOR AND SPATIAL DIMENSION. The correction factor _D is the ratio of the mean error over the square root of the mean square error. Its theoretical value depends only on the dimension D of the work space, assuming Gaussian error distribution (APPENDIX B). As a consequence, its empirical value could be used to estimate the dimension D. For reconstruction errors obtained for all time steps by the Bayesian two-step method, the empirical value was 0.8705 for animal 1 and 0.8855 for animal 2, surprisingly close to the value of the 2-D correction factor ₂ = 0.8862 but quite different from ₁ and ₃ (Table B1). However, the actual error distribution differed from the theoretical probability density ₂(r) (Table B1 and Fig. 16) because it peaked at a smaller value than the standard deviation (r = s) and had a longer tail for large errors. This might be related to the fact that the rats run in 2-D space but within the confinement of a track so that the work space was neither purely two dimensional nor purely one dimensional.

COUNTERINTUITIVE EFFECT OF PLACE FIELD SIZE. Imagine we have a fixed number of place cells all having place fields within a given region. Suppose we can change the size of the place fields without altering the positions of the place field centers and the peak firing rates. How would this change affect the accuracy of spatial information encoded in the population? The same question in a different form is how the minimal achievable reconstruction error depends on the place field size. This question is related to how the accuracy of spatial representation is affected by the broader place fields found in ventral hippocampus of rats (Jung et al. 1994) and in the hippocampus of the genetically engineered mouse (see further). The size of place fields also tends to increase during repetitive running on a track (Mehta and McNaughton 1997).

When the peak firing rate is fixed, changing the size of place fields should have no effect on the minimal reconstruction error or the accuracy of spatial representation. This is because Eq. 56 does not contain the width  of the place fields. Here we have assumed that the place fields are of Gaussian shape, the spikes are Poisson, and different cells are independent.

The explanation is that larger place fields imply more overlaps and therefore higher overall average firing rate of the whole population; this will reduce reconstruction error. This effect turns out to precisely compensate for the less accurate spatial information provided by a single enlarged place field. A confirmation by simulation is shown in Fig. 10. The fact that larger place field sizes imply higher mean firing rates is clearly seen in Eq. 55. This effect precisely cancels out the effect of the tuning width  in Eq. 53 so that the width  is absent in the equivalent Eq. 56.

A related issue is the recent reports of increased place field size for CA1 cells in the hippocampus of genetically altered mice with impairment of long-term potentiation (LTP) restricted to CA1 (McHugh et al. 1996; Rotenberg et al. 1996). A decease of the peak firing rate of place cells was reported explicitly by Rotenberg, Mayford, Hawkins, Kandel, and Muller (1996). According to Eq. 55 such a decrease also was implied in the observation that the mean firing rate of the genetically altered mice were the same as the normal (McHugh et al. 1996). The increase of reconstruction error for these mice (McHugh et al. 1996) may be attributed to a decrease in peak firing rate rather than to an increase in place field size per se because Eq. 56 suggests that reconstruction error should be independent of place field size. Equivalently, if the mean firing rate is fixed, it is equally valid to apply Eq. 53, which suggests that in this case the error should be larger for larger place fields. Instability of place fields is another possible source of reconstruction error (cf. Rotenberg et al. 1996). A more detailed analysis of the data would be required for quantitative comparison.

The accuracy of spatial representation can be improved by increasing the number of cells or their firing rates but not by sharpening the place fields. For example, after repetitive running on a track, both the size of the place fields and the average firing rate tended to increase (Mehta and McNaughton 1997). The increased firing rate would imply better accuracy of spatial representation, regardless of the place field size.

General biological discussion

The minimal achievable reconstruction considered in the preceding sections is actually a general measure of how accurately a physical variable is encoded by a population of neurons, regardless of the decoding method used. As a consequence, the theoretical limit is valid not only for the computational problem of reconstruction but also for all neurophysiological and behavioral measures of acuity or accuracy in the sense of mean square error, as long as the information is derived solely from the same neuronal population.

We have derived explicit formulas for cells with Gaussian and cosine tuning functions because of their common usage, although the methods can be adapted readily for more general situations. From these formulas, we know precisely how the reconstruction accuracy depends on the parameters of individual cells and the cell number density. For example, suppose we have spatial region of 1 × 1 m. How many place cells must be used to achieve a spatial acuity of d = 1 cm everywhere in this region? Using Eq. 56, we have

(72)

Assuming a mean peak firing rate f_max ~ 15 Hz and a biologically plausible time window  ~ 200 ms, we get N ~ 10³. That is, a minimum of ~1,000 cells would be required. As another example, using up all the place cells in a rat, say, N ~ 10⁵ (Amaral and Witter 1995), what is the maximum area that can be covered under the acuity of 1 cm? Using the same parameters considered earlier, we find the total area is A ~ 10⁶ cm², the same as a 10 × 10 m square.

The insensitivity of place field size on reconstruction accuracy is a peculiarity of two dimensions. As shown in Fig. 10, more information is encoded by a sharper tuning curve in 1-D space. Examples of 1-D tuning curves include the orientation tuning in visual cortex and directional tuning of head-direction cells (Taube et al. 1990). In contrast, in3-D space or space of higher dimensions, a sharper tuning width implies worse information encoding. So the 2-D example of place cell is a critical case where the tuning width does not matter. To increase representation accuracy, one should use neurons with broader tuning functions to represent a variable higher than two dimensions, assuming that the peak firing rate is fixed. However, broader tuning width implies higher energy consumption because more neurons are activated at any given time. See Hinton, McClelland, and Rumelhart (1986), Baldi and Heiligenberg (1988), Zohary (1992), Gerstner, Kempter, van Hemmen, and Wagner (1996) for related discussions on the consequences of broad tuning.

	BIOLOGICAL PLAUSIBILITY OF RECONSTRUCTION METHODS

Unified formulation

The main result in this section is that all the reconstruction methods discussed in this paper can be implemented as linear feedforward neural networks. This result is in sharp contrast with the general belief that reconstruction methods such as the Bayesian method or the template matching method are merely mathematical techniques without any biological relevance.

The reconstructed position in all reconstruction methods considered above can be expressed as

(73)

and function (x) can be written in the form

(74)

where as before n_i is the number of spikes that cell i fired in a certain time window, _i(x) is a basis function associated with cell i. The function A(x) is an additive bias, which is independent of the activity of the cells and is reminiscent of a regularization term (Poggio et al. 1985; Tikhonov and Arsenin 1977). See also Zemel, Dayan, and Pouget (1997) for related discussion. The bias is needed only for implementing Bayesian methods. For all basis function methods, there is no bias term; that is, A(x)  0. See Table 2 for a summary.

Here the Bayesian methods are implemented by taking the logarithm of the posterior probability in Eqs. 36 and 41. Because logarithm is monotonic and does not change the peak position, maximizing the product f_i(x)^n_i is equivalent to maximizing the sum n_ii(x), where

(75)

To avoid taking the logarithm of zero, an arbitrary small positive number may be added to the tuning function f_i. The other factors in the Bayesian Eq. 36 becomes the additive bias

(76)

which does not depend on activity n_i. An intuitive interpretation of the bias is that it favors positions that animal visits frequently, that is, places with high P(x), and at the same time disfavors places that are over-represented by a large number of cells or cells with excessively high firing rates. The bias term becomes a constant and can be ignored if the animal is equally likely to visit all positions and the average firing rates of all cells are uniform across space.

The difference between addition and multiplication of basis functions is illustrated in Fig. 13, although taking logarithm can transform a product into a sum. For the direct basis method, the computation is the sum

(77)

In Bayesian method, the crucial part of the main Eq. 36 is the product

(78)

The product is like a logical "and" operation and the final peaks correspond to those positions where all contributing template functions in the product are not close to zero. By contrast, the sum behaves like a logical "or" operation and all peaks from contributing template functions are preserved. This difference explains why output distributions of basis function methods are usually much broader and more likely to contain multiple peaks than that of the Bayesian methods (Fig. 2).

View larger version (12K): [in this window] [in a new window]   FIG. 13. An example of addition and multiplication of 2 basis functions. Addition yields a distribution with 2 peaks, yet multiplication yields only 1 peak. Key operation of Bayesian reconstruction method is multiplication and its output distribution is often sharper and less likely to contain multiple peaks than output distribution of basis function method, whose key operation is addition. See Fig. 2 for an actual example of comparison.

Biological implications

The unified formulation of reconstruction methods in the preceding section can be implemented readily by a biologically plausible neural network structure. The first computational step in Eq. 74 is linear combination of fixed basis functions. Here the coefficients are the numbers of spikes, which are proportional to the instantaneous firing rates within the given time window, which is plausibly a fraction of a second. The linear combination only requires a feedforward network in which the basis functions are implemented by the connection weights (Fig. 14). More precisely, the sum n_ii(x) is computed as

(79)

View larger version (10K): [in this window] [in a new window]   FIG. 14. A unified network implementation for various reconstruction methods. First layer contains N cells tuned arbitrarily to an variable x. Cells in 2nd layer represent value of x by their locations in layer. Activity distribution in 2nd layer represents how likely it is for each possible value of x to be true value.

Indeed, linear combination of basis functions is biologically plausible computational strategy that has been used by many authors, such as in 3-D object representations (Poggio and Girosi 1990; Ullman and Basri 1991), parietal spatial representation (Pouget and Segnowski 1994, 1997; Zipser and Andersen 1988), modeling place cells (O'Keefe and Burgess 1996; Redish and Touretzky 1996; Zipser 1985), and/or as a general strategy for neural representation (Anderson and Van Essen 1994; Poggio 1990; Zemel et al. 1997).

The maximization in the second computational step in Eq. 73 can be implemented by any circuits that can approximate a winner-take-all operation on the activity distribution in the second layer after the first feedforward step. An exact winner-take-all implementation might not be needed in the biological system because there are different ways of reading the distributed information after the first step is done. In general, the full activity distribution after the first feedforward step contains more information about the probability distribution of the position and it may be useful to maintain this information for the purpose of further computations that need to take into account the variance of the estimate as well as the mean of the position (Nowlan and Sejnowski 1995).

The bias term for the one-step Bayesian method can be implemented simply as constant inputs to the second layer. To implement the additional bias for the two-step Bayesian method, a type of short-term memory is needed for the winner-take-all operation. The quadratic bias in Table 2 comes from the Gaussian distribution in Eq. 43, but the precise functional form is not critical. In network implementation, any facilitation mechanism that favors units close to the previously selected winner unit and prohibits units far away should suffice.

Therefore, it is feasible for a biological system to implement all the reconstruction methods considered in this paper. Moreover, a simple Hebb rule should be enough for learning the weights used in the Bayesian method and the direct basis method because the desired template function for each cell depends only the tuning function of the same cell. However, the reciprocal basis method may need a nonlocal learning rule.

A related observation is that attractor dynamics in a recurrent network may perform operations similar to winner-take-all, noise cleaning or maximum likelihood estimation (Pouget and Zhang 1997). This approach has some biological plausibility because attractor dynamics might actually be used by the brain, for example, in the hippocampal system (McNaughton et al. 1996; Samsonovich and McNaughton 1997; Tsodyks and Sejnowski 1995a; Zhang 1996), the primary visual cortex (Ben-Yishai et al. 1995; Somers et al. 1995; Tsodyks and Sejnowski 1995b), and motor cortex (Lukashin et al. 1996).

One may ask which reconstruction algorithm is actually used by the biological system, for example, for spatial navigation with the hippocampus place cells. There are a few general predictable effects. According to the unified network implementation in Eq. 74, different methods are distinguished mainly by the template functions that are implemented as synaptic weights. For example, the direct basis method would imply a weight pattern directly proportional to the mean firing rate given by the tuning function, whereas the Bayesian method would imply a weight pattern proportional to the mean firing rate on a logarithm scale. Whether or not the brain uses the most efficient construction method might also be tested by spatial acuity at neural or behavioral levels. For our dataset, such information is not available. What we have shown here is that the biological system does have the computational resources to implement the most efficient algorithms we considered, which can reach the best possible accuracy defined by Cramér-Rao lower-bound under idealized situations.

	CONCLUSIONS

The goal of this paper has been in part to compare different reconstruction methods, to improve the accuracy of these methods, and to assess their performance against the theoretical lower bound on reconstruction accuracy for all possible methods. For trajectory reconstruction based on spike trains from simultaneously recorded hippocampal place cells, probabilistic methods based on Bayesian formula were exceptionally accurate, especially when information about the reconstructed position from the previous time step was used to discourage discontinuities in the trajectory. In our best reconstruction with an average of up to ~30 spikes within the time window, the mean errors were in the range of the 5 cm intrinsic error of the position tracking system.

This reconstruction study has also revealed some interesting properties of place cells, such as erratic jumps of reconstruction errors and their correlation with movement and the local field potentials, the slow systematic drift of reconstruction error, as well as a possible shift of the optimal center of the reconstruction time window towards the past during continuous periodic running.

There are two major implications for our consideration of the theoretical lower bound on reconstruction accuracy derived by Fisher information and Cramér-Rao inequality. First, this lower bound determines precisely how accurately a physical variable is encoded in the neuronal population in the sense of mean square error, regardless of which method is used for reading out the information. That is, this bound quantifies the amount of information encoded in the neuronal population regardless of decoding method. As a consequence, the best achievable accuracy is valid not only for the reconstruction problem, but also for all possible neural or behavioral measures of acuity or accuracy, as long as the information is derived solely from the neuronal population in question. Second, we have identified that the Bayesian reconstruction methods can reach the theoretical lower bound when the conditions of Poisson spike statistics and independence of different cells are satisfied. In such situations, optimal reconstruction has been achieved. For our place cell dataset, the best reconstruction errors were above the estimated theoretical limit by perhaps a factor of two or so. The remaining discrepancy may be caused by various approximations involved.

One counterintuitive result emerging from our theoretical consideration of accuracy is that making the size of all place fields smaller (or larger) without changing their peak firing rate would have no effect on the accuracy of coding the animal's position in the whole population. However for1-D work space, a narrower tuning curve gives more accurate coding. By contrast, for three and higher dimensions, a broader tuning function gives more accurate coding. Thus two dimensions are a critical case in which the coding accuracy is independent of the width of the tuning function.

Despite popular belief to the contrary, we have shown that all reconstruction methods considered in this paper including template matching and the Bayesian method can be implemented as simple feedforward neural networks. This result suggests that biological systems are capable of implementing all of these reconstruction methods. The reduction of the nonlinear Bayesian method to the linear basis function framework relies on the assumptions of Poisson spiking statistics and independence of different cells. For the place cell data, they are reasonable approximations, although they are not exactly true.

Whenever neuronal activity is correlated with a measurable physical variable, reconstruction of this variable from population activity is a relevant problem both as a research tool and as a hint to how the brain might solve the same problem. There are no intrinsic constraints on the type of physical variables that can be reconstructed, or on the type of tuning functions that the cells can encode. For example, instead of the position in a continuous 2-D space, the variable could as well be discrete and disjoint, which might be more suitable for representing distinct classes of objects or categories.

	ACKNOWLEDGEMENTS

  We thank H. Kudrimoti and W. E. Skaggs for valuable help in data processing and for many discussions, J. Wang for cluster cutting, M. S. Lewicki for discussions on Bayesian methods, A. Pouget for discussions on population coding and Fisher information, A. B. Schwartz and R. U. Muller for helpful suggestions on the manuscript, and two anonymous reviewers for useful comments. Data acquisition by J. L. Gerrard and C. A. Barnes.

  This work was supported by Howard Hughes Medical Institute and National Institute of Neurological Disorders and Stroke Grant NS-20331.

	APPENDIX A: RECIPROCAL BASIS WITH MODIFICATIONS

In this section we modify the reciprocal basis by taking into account a nonuniform probability distribution for visiting different positions as well as the randomness in firing rate. As with the original reciprocal basis, the modified basis functions are optimal with respect to mean square error. However the reciprocal relationship is lost. For our dataset, the performances of these modified methods were sometimes better than the original reciprocal basis, as shown in Table 1.

In general, consider a basis function method

(A1)

where N is the total number of cells, n = (n₁, . . . , n_N) are the numbers of spikes fired by these cells within the time window  and (g₁, . . . , g_N) are the basic functions to be determined. Here  is included in the formula for the convenience of scaling; it does not change the result. The discrete version of the basis function g_i(x) is the column vector g_i, which is defined by first digitizing the position x and then concatenating the pixels. The vectors are collected into a single matrix G = [g₁, . . . , g_N]. As shown in Eq. 15, the simple reciprocal basis is related to the firing rate matrix F by

(A2)

where the second equality is valid only when the matrix inverse is nonsingular.

The reciprocal basis functions obtained by Eq. A2 are appropriate only if the animal visits each spatial position equally frequently. When the visiting probability is nonuniform, there are several ways to modify the matrix G or the basis functions. First the simple Eq. 21 is the same as the prescription

(A3)

This is the formula used in our systematic comparison of different reconstruction methods in Fig. 5. Alternatively, we may choose

(A4)

A third possible formulation is

(A5)

Here three diagonal matrices have been used, which are defined as

(A6)

(A7)

and

(A8)

where K is the total number of position pixels, P_k is the probability for visiting position k, and _i is the average firing rate of cell i.

These formulas are related. In Eq. A5, the D term is related to firing rate variability (see below). For large time window , the D term vanishes, and Eq. A5 is reduced to Eq. A4. When the visiting probability P_k is the same for all positions, Eq. A4 is reduced to the original reciprocal basis Eq. A2.

Equation A5 has been derived by minimizing the error function

(A9)

(A10)

When firing rate variability is ignored, the equation

(A11)

holds exactly. This is a reasonable approximation when the firing rates are high or the time window is large. Dropping out the term P_n|k in Eq. A9, we only need to minimize

(A12)

The result is Eq. A4.

	APPENDIX B: MINIMAL RECONSTRUCTION ERROR BASED ON FISHER INFORMATION

Basic approach

We provide an intuitive argument for using Fisher information and then sketch our method for estimating a lower bound for reconstruction error using the Cramér-Rao inequality.

As an intuitive example, consider a 1-D situation with an idealized place cell whose firing rate f(x) depends only on position x (Fig. B1). Suppose the cell fired n spikes within the time window . Because of random variation in the exact number of spikes, the formula

(B1)

is true only on average but not for each individual trial, When we use this single cell to estimate the position x from the actual firing rate n/ and the shape of the tuning curve f(x), the fluctuation of the number of spikes would result in an error x that obeys

(B2)

Thus the mean squared error is

(B3)

In the last step, Poisson spiking statistics was assumed; that is, the variance n² is equal to the mean n = f(x). Here the reconstruction error is expressed in terms of the tuning function f(x) and its slope f(x). It makes sense that a larger slope implies smaller error because a small change of position would imply a large and noticeable change in the firing rate (Fig. B1).

View larger version (12K): [in this window] [in a new window]   FIG. B1. Intuitive argument for Fisher information in estimating position x from instantaneous firing rate and a known tuning curve f(x). Fluctuation in firing rate leads to an error x in estimated position, which decreases with slope of tuning curve.

An identical result is obtained by using the formal definition of Fisher information (Cover and Thomas 1991; Kay 1993; Rao 1965; Scharf 1991; Seung and Sompolinsky 1993; Zacks 1971):

(B4)

where the last equality follows from straightforward calculation under the assumption that the probability P(n|x) of firing n spikes at position x has Poisson distribution with the mean firing rate f(x).

The squared error for any unbiased estimator of position x cannot exceed the Cramér-Rao bound, namely

(B5)

whose validity originates from the universal Cauchy-Schwartz inequality. This result agrees with the intuitive Eq. B3.

When more than one cell is used in reconstruction, as long as they are independent, the total Fisher information is the sum

(B6)

This leads to an effective composition rule of minimal errors:

(B7)

where ² is the minimum achievable value of x² by using all the cells from 1 to N, and ²_i is defined as the same error using only cell i.

So far, we have only considered a 1-D problem. Extension to 2-D estimates is simple if we assume that the error statistics in x and y directions are the same. Then the squared error of reconstruction is

(B8)

where J is the Fisher information for variable x alone. The minimal reconstruction error is estimated as ₂. The final result is given by Eq. 56 or Eq. 53.

Derivation of error formula in arbitrary dimensions Following the same approach considered in the preceding section, we list the keys steps in deriving the general formula for minimal reconstruction error in arbitrary spatial dimension D. The mean firing rate of each cell is assumed to obey a Gaussian tuning function in a D-dimensional space

(B9)

(B10)

The Fisher information with respect to a single dimension, x₁ for example, is

(B11)

where the first equality is a definition, the second equality follows from the Poisson statistics Eq. B10 and the last equality follows from the Gaussian tuning Eq. B9.

First consider the simple case where all cells have identical peak firing rate and tuning width. Now the total Fisher information with respect to x₁ for all N cells should be

(B12)

where c_i is the center of the Gaussian tuning function for cell i, and  is the number of cells per unit volume in D-dimensional space (see Eq. 62). The sum is replaced by the integral under the assumption that the number of cells N is large and the centers c₁, . . . , c_N are uniformly distributed in space. Because of the symmetry with respect to different dimensions, we can use the same argument in Eq. B8 to obtain the Fisher information for all N cells with all dimensions included

(B13)

If different cells have different peak firing rate f_max and tuning width , one can easily see that Eq. B13 should be modified to

(B14)

where the average   is with respect to the population of cells with different peak firing rate and tuning width. Assuming independence of the tuning parameters, we have

(B15)

Finally

(B16)

where the correction factor F_D is given by Eq. B26. The result is Eq. 58. In addition, it is easy to derive Eq. 61, which immediately leads to the equivalent Eq. 57.

Estimating tuning width To evaluate the Fisher information, the width  in the Gaussian tuning model needs to be estimated. Given the spatial distribution of the average firing rate f(x, y) measured experimentally for a single cell, we want to estimate the tuning width  assuming a Gaussian model

(B17)

There are several methods for estimating . The first method uses entropy or Shannon information

(B18)

This method requires an estimate of the peak rate f_max, which is sensitive to noise and smoothing.

The second method uses the variance

(B19)

In the above, all averages   should be interpreted as follows. Suppose each unit square in the grid has the area h² and f(i, j) is the value of average firing rate at grid point (i, j), then f = h²f(i, j) and so on. Equations B18 and B19 are closely related to the information content measure and sparsity measure considered by Skaggs et al. 1996. One difference is that here the probability for spatial occupancy should not be included.

The third method exploits the linear relationship between the area A of a place field above a given cutoff firing rate f_cutoff and the logarithm of f_cutoff (Muller et al. 1987). This linearity is consistent with the Gaussian model in Eq. B17, which gives

(B20)

Because A depends linearly on ln f_cutoff, both  and f_max can be obtained by linear regression.

In our dataset, the movement of the animal was restricted on a narrow track. Simply assuming that each observed place field is a circular symmetric 2-D Gaussian cut by the track, we can estimate the width  of the uncut Gaussian by a modified variance method. First, concatenate the pixels in the 2-D distribution f(i, j) to obtain a 1-D distribution f(i) and then assume a 1-D Gaussian distribution to obtain the 1-D tuning width

(B21)

Here f = hf(i) and so on, with h being the interval in the1-D grid. The 2-D tuning width  should equal approximately ₁ divided by the number of grid points spanned by the width of the track. This is because when a 2-D Gaussian is cut by an arbitrary straight line, the resultant 1-D distribution is a Gaussian of the same tuning width. The above procedure is valid if different 1-D slices cut in parallel with the track have approximately the same peak rate.

Distribution of reconstruction error and correction factor A correction factor is needed because the Cramér-Rao bound only gives an estimate of the variance of unbiased reconstruction, that is, the average square of the reconstruction error r². Its square root differs slightly from the direct average of reconstruction error r. The task of this section is to determine the ratio _D = r/ analytically.

For reconstruction in D-dimensional space, it is natural to assume that the distribution of the error vector obeys a D-dimensional Gaussian distribution

(B22)

(B23)

The probability density _D(r) for reconstruction error can thus be obtained from Eqs. B22 and B23:

(B24)

(B25)

For example, when m = 2 the formula has a particularly simple form r² = Ds², and (r/s)² obeys the chi square distribution with D degrees of freedom.

The correction factor, defined as r/, is thus

(B26)

Note that _D depends only on the spatial dimensionD. The width parameter s for the error distribution has disappeared. Although _D < 1 for all dimension D, it approaches 1 monotonically as D  .

The results for 1-, 2-, and 3-D spaces are listed in Table B1. The theoretical curves can accurately describe the empirical error distribution of Bayesian reconstruction on synthetic data (Fig. B2). This result is consistent with the assumption of Gaussian error distribution in the original D-dimensional space. In fact, for a locally uniform prior probability distribution, Bayesian reconstruction method is equivalent to maximum likelihood estimation, the error distribution of which is known to be asymptotically Gaussian in the limit of large sample (Cramér 1946).

View larger version (17K): [in this window] [in a new window]   FIG. B2. Distributions of reconstruction errors in 1-, 2-, and 3-D spaces. Continuous curves are theoretical distributions in Table B1, which agree nicely with histograms obtained by simulations of Bayesian reconstructions repeated for 5,000-10,000 trials, assuming Gaussian tuning function and Poisson spike statistics. Tuning width was  = 3. All other parameters were identical to those in Fig. 10.

	APPENDIX C: OPTIMALITY OF BAYESIAN RECONSTRUCTION

As shown in Figs. 10 and 12, Bayesian reconstruction in Eq. 36 can achieve the best possible accuracy in the sense of mean square error defined by Fisher information and the Cramér-Rao bound, provided that the conditions of Poisson spike statistics and independence of different cells are satisfied. In such a situation, the optimality of Bayesian reconstruction method can be verified directly as follows. This analysis can provide an explicit estimate of the reconstruction error in each individual trial. In general, for a uniform prior probability distribution, Bayesian reconstrution is equivalent to maximum likelihood estimation, which can reach the Cramér-Rao bound in the limit of large sample (Cramér 1946).

Consider the Fisher information for a 1-D reconstruction under a uniform prior probability:

(C1)

where f_i(x) is the tuning curve for cell i and  is the time window (see APPENDIX B). Let x be the reconstruction error of the Bayesian method, which seeks the position x that maximizes

(C2)

(C3)

(C4)

when a large number of cells is used. That is, when the Cramér-Rao lower bound is achieved.

Consider how fluctuation in the instantaneous firing rate n_i/ affects the estimated x because without the fluctuation the reconstruction would always be perfect. By definition, the average of n_i/ is f_i(x). Calling the difference _i, we write

(C5)

Replace x by x + x in Eq. C3 and then make Taylor expansion, keeping only the first-order term of x . After straight forward algebra we find

(C6)

This shows how the estimated position depends on the firing rate fluctuation in an individual trial. The above approximation procedure is valid if a large number of cells is used so that the error is small. Because by definition _i = 0, it follows from Eq. C6 that the estimated position must be unbiased, i.e., x = 0. Squaring Eq. C6 and averaging yield

(C7)

where all cross terms _ij with i  j have vanished because of the independence of different cells. By using the definition Eq. C1 and Poisson spike statistics ²_i = f_i(x)/, we obtain Eq. C4 from Eq. C7 as desired.

	APPENDIX D: ERROR OF POPULATION VECTOR UNDER COSINE TUNING

2-D work space

Consider N cells with cosine tuning functions and preferred rections scattered randomly around the circle with a uniform distribution. Let the spikes of different cells occur independently of each other and obey a Poisson distribution. The task is to estimate the average angular error of the population vector.

Note that by average over trials, denoted by  _T, we mean that the preferred directions of all cells are fixed and the average is performed over repetitive trials where the only randomness is the Poisson statistics of the spikes; by average over configurations, denoted by  _C, we mean averaging over random choices of the preferred directions. The error of population vector comes from both the randomness of the spikes and the randomness of the preferred directions.

Because of the rotation symmetry, we need to consider only a single test direction along the x-axis, u = (1, 0). The average firing rate of cell i is

(D1)

where A and B are constants (cf. Eq. 64) and vector u_i = (u_ix, u_iy) is of unit length and represents the preferred direction of cell i. The population vector is defined as

(D2)

with the coefficients given by

(D3)

where n_i is the number of spikes for cell i within the time window  and the background firing rate B is subtracted to improve reconstruction accuracy (Georgopoulus et al. 1988). The x component of the population vector is estimated as follows

(D4)

in which the two approximation steps are valid if the number of cells N is large. In the above derivation, we have evaluated several averages

(D5)

follows directly from the definition Eq. D3 and the fact that n_iT = f_i; the equations

(D6)

are a consequence of the uniform distribution of preferred direction around the circle.

Because the test direction is assumed to be along the x-axis, population vector U is also close to the x-axis. Thus |U|  U_x and the angular error of the population vector can be approximately as

(D7)

where U_x is a constant given by Eq. D4 and U_y =  c_iu_iy with c_i given by Eq. D3. It follows immediately that

(D8)

After expanding U²_yT, evaluate each term using

(D9)

which follows from Poisson spike statistics and the independence of different cells. Next, substitute n_iT = f_i = (Au_ix + B) into the expansion and note that all terms containing odd powers such as u_ixu²_iy and u_ixu_jxu_iyu_jy will vanish after averaging over configurations. Thus we need only keep terms with even powers and obtain

(D10)

which can be evaluated by applying Eq. D6 again. The final result of the average angular error is ₁ = ₁ with Q₁ given by Eq. 69.

3-D work space We need to consider only a single test direction along the x-axis: u = (1, 0, 0). The method is very similar to that in the case of 2-D work space, except that now

(D11)

In addition, the error now occurs in both y and z directions so that

(D12)

The final result of the average angular error is ₂ = ₂ with Q₂ given by Eq. 70.

	FOOTNOTES

¹   For comparison, we tested square root law on data taken from Fig. 3 in paper by Georgopoulos, Kettner and Schwartz (1988), where error of population vector was computed by using up to 475 motor cortical cells. Errors were given in terms of 95% confidence cone, which was expected to be proportional to mean angular error. Square root law held quite well: on log-log scale, a linear fit yielded a slope of 0.519 with a correlation coefficient of 0.9984, as compared with theoretical slope of 1/2.

  Address reprint requests to K. Zhang.

  Received 27 May 1997; accepted in final form 25 August 1997.

	REFERENCES

Abstract Introduction References

• ABBOTT, L. F. Decoding neuronal firing and modelling neural networks. Q. Rev. Biophys. 27: 291-331, 1994.[Medline]
• ABBOTT, L. F., BLUM, K. I. Functional significance of long-term potentiation for sequence learning and prediction. Cereb. Cortex 6: 406-416, 1996.[Abstract/Free Full Text]
• AMARAL, D. G., WITTER, M. P. Hippocampal formation. In: The Rat Nervous System, edited by and G. Paxinos , 2nd ed.. San Diego, CA: Academic, 1995, p. 443-493
• ANDERSON, C. H. AND VAN ESSEN, D. C. Neurobiological computational systems. In: Computational Intelligence: Imitating Life, edited by J. M. Zurada, R. J. Marks II, and C. J. Robinson. New York: IEEE Press, 1994, p. 213-222.
• BALDI, P., HEILIGENBERG, W. How sensory maps could enhance resolution through ordered arrangements of broadly tuned receivers. Biol. Cybern. 59: 313-318, 1988.[Medline]
• BARNES, C. A., SUSTER, M. S., SHEN, J., MCNAUGHTON, B. L. Multistability of cognitive maps in the hippocampus of old rats. Nature 388: 272-275, 1997.[Medline]
• BEN-YISHAI, R., BAR-OR, R. L., SOMPOLINSKY, H. Theory of orientation tuning in visual cortex. Proc. Natl. Acad. Sci. USA 92: 3844-3848, 1995.[Abstract/Free Full Text]
• BIALEK, W., RIEKE, F., DE RUYTER VAN STEVENINCK, R. R., WARLAND, D. Reading a neural code. Science 252: 1854-1857, 1991.[Abstract/Free Full Text]
• BLAIR, H. T., LIPSCOMB, B. W., SHARP, P. E. Anticipatory time intervals of head-direction cells in the anterior thalamus of the rat: implications for path integration in the head-direction circuit. J. Neurophysiol. 78: 145-159, 1997.[Abstract/Free Full Text]
• BLAIR, H. T., SHARP, P. E. Anticipatory head-direction cells in anterior thalamus: evidence for a thalamocortical circuit that integrates angular head velocity to compute head direction. J. Neurosci. 15: 6260-6270, 1995.[Abstract]
• BLUM, K. I, ABBOTT, L. F. A model of spatial map formation in the hippocampus of the rat. Neural Comput. 8: 85-93, 1996.[Medline]
• BROWN, E. N., FRANK, L., WILSON, M. Statistical approaches to place field estimation and neuronal population decoding. Soc. Neurosci. Abstr. 22: 910, 1996.
• BURGESS, N., RECCE, M., O'KEEFE, J. A model of hippocampal function. Neural Networks 7: 1065-1081, 1994.
• BUZSÁKI, G. Hippocampal sharp waves: their origin and significance. Brain Res. 398: 242-252, 1986.[Medline]
• BUZSÁKI, G. Two-stage model of memory trace formation: a role for "noisy" brain states. Neuroscience 31: 551-570, 1989.[Medline]
• COVER, T. M., THOMAS, J. A. In: Elements of Information Theory., . New York: Wiley, 1991
• CRAMÉR, H. In: Mathematical Methods of Statistics., . Princeton, NJ: Princeton Univ. Press, 1946
• EICHENBAUM, H., WIENER, S. I., SHAPIRO, M. L., COHEN, N. J. The organization of spatial coding in the hippocampus: a study of neural ensemble activity. J. Neurosci. 9: 2764-2775, 1989.[Abstract]
• FENTON, A. A. AND MULLER, R. U. Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proc. Natl. Acad. Sci. USA. In press.
• FÖLDIÁK, P. The "ideal humunculus": statistical inference from neural population responses. In: Computation and Neural Systems 1992, edited by F. Eeckman, and J. Bower . Norwell, MA: Kluwer, 1993
• GEORGOPOULOS, A. P., KETTNER, R. E., SCHWARTZ, A. B. Primate motor cortex and free arm movements to visual targets in three-dimensional space. II. Coding of the direction of movement by a neuronal population. J. Neurosci. 8: 2928-2937, 1988.[Abstract]
• GEORGOPOULOS, A. P., LURITO, J. T., PETRIDES, M., SCHWARTZ, A. B., MASSEY, J. T. Mental rotation of the neuronal population vector. Science 243: 234-236, 1989.[Abstract/Free Full Text]
• GEORGOPOULOS, A. P., SCHWARTZ, A. B., KETTNER, R. E. Neuronal population coding of movement direction. Science 233: 1416-1419, 1986.[Abstract/Free Full Text]
• GERRARD, J. L., GINZBURG, I., BARNES, C. A., SEJNOWSKI, T. J., MCNAUGHTON, B. L. Accuracy of hippocampal place fields in predicting location is enhanced in the presence of theta rhythm. Soc. Neurosci. Abstr. 21: 942, 1995.
• GERSTNER, W., KEMPTER, R., VAN HEMMEN, J. L., WAGNER, H. A neuronal learning rule for sub-millisecond temporal coding. Nature 383: 76-78, 1996.[Medline]
• GOLUB, G. H., VAN LOAN, C. F. In: Matrix Computations, , 3rd ed.. Baltimore, MD: Johns Hopkins Univ. Press, 1996
• HINTON, G. E., MCCLELLAND, J. L., RUMELHART, D. E. Distributed representations. In: Parallel Distributed Processing, edited by D. E. Rumelhart, and J. L. McClelland . Cambridge, MA: MIT Press, 1986, vol. 1, p. 77-109
• JENSEN, O., LISMAN, J. E. Hippocampal CA3 region predicts memory sequences: accounting for the phase precession of place cells. Learning and Memory 3: 279-287, 1996. [Abstract/Free Full Text]
• JUNG, M. W., WIENER, S. I., MCNAUGHTON, B. L. Comparison of spatial firing characteristics of units in dorsal and ventral hippocampus of the rat. J. Neurosci. 14: 7347-7356, 1994.[Abstract]
• KAY, S. M. In: Fundamentals of Statistical Signal Processing: Estimation Theory., . Upper Saddle River, NJ: Prentice Hall, 1993
• KUBIE, J. L., FOX, S. E., MULLER, R. U. Variations in place cell firing with the state of the hippocampal EEG. Soc. Neurosci. Abstr. 10: 599, 1984.
• LEHKY, S. R., SEJNOWSKI, T. J. Neural model of stereoacuity and depth interpolation based on a distributed representation of stereo disparity. J. Neurosci. 10: 2281-2299, 1990.[Abstract]
• LUKASHIN, A. V., AMIRIKIAN, B. R., MOZHAEV, V. L., WILCOX, G. L., GEORGOPOULOS, A. P. Modeling motor cortical operations by an attractor network of stochastic neurons. Biol. Cybern. 74: 255-261, 1996.[Medline]
• MCHUGH, T. J., BLUM, K. I., TSIEN, J. Z., TONEGAWA, S., WILSON, M. A. Impaired hippocampal representation of space in CA1-specific NMDAR1 knockout mice. Cell 87: 1339-1349, 1996.[Medline]
• MCNAUGHTON, B. L., BARNES, C. A., GERRARD, J. L., GOTHARD, K., JUNG, M. W., KNIERIM, J. J., KUDRIMOTI, H., QIN, Y., SKAGGS, W. E., SUSTER, M., WEAVER, K. L. Deciphering the hippocampal polyglot: the hippocampus as a path integrating system. J. Exp. Biol. 199: 173-185, 1996.[Abstract]
• MCNAUGHTON, B. L., BARNES, C. A., O'KEEFE, J. The contributions of position, direction, and velocity to single unit activity in the hippocampus of freely-moving rats. Exp. Brain Res. 52: 41-49, 1983.[Medline]
• MEHTA, M. R., BARNES, C. A., MCNAUGHTON, B. L. Experience-dependent, asymmetric expansion of hippocampal place fields: Proc. Natl. Acad. Sci. USA 94: 8918-8921, 1997.[Abstract/Free Full Text]
• MULLER, R. U., BOSTOCK, E., TAUBE, J. S., KUBIE, J. L. On the directional firing properties of hippocampal place cells. J. Neurosci. 14: 7235-7251, 1994.[Abstract]
• MULLER, R. U., KUBIE, J. L. The firing of hippocampal place cells predicts the future position of freely moving rats. J. Neurosci. 9: 4101-4110, 1989.[Abstract]
• MULLER, R. U., KUBIE, J. L., RANCK, J. B. JR Spatial firing patterns of hippocampal complex-spike cells in a fixed environment. J. Neurosci. 7: 1935-1950, 1987.[Abstract]
• NOWLAN, S., SEJNOWSKI, T. J. A selection model for motion processing in area MT of primates. J. Neurosci. 15: 1195-1214, 1995.[Abstract]
• O'KEEFE, J., BURGESS, N. Geometric determinants of the place fields of hippocampal neurons. Nature 381: 425-428, 1996.[Medline]
• O'KEEFE, J., DOSTROVSKY, J. The hippocampus as a spatial map: preliminary evidence from unit activity in the freely-moving rat. Brain Res. 34: 171-175, 1971.[Medline]
• O'KEEFE, J., RECCE, M. L. Phase relationship between hippocampal place units and the EEG theta rhythm. Hippocampus 3: 317-330, 1993.[Medline]
• OLTON, D. S., BRANCH, M., BEST, P. J. Spatial correlates of hippocampal unit activity. Exp. Neurol. 58: 387-409, 1978.[Medline]
• ONO, T., NAKAMURA, K., NISHIJO, H., EIFUKU, S. Monkey hippocampal neurons related to spatial and nonspatial functions. J. Neurophysiol. 70: 1516-1529, 1993.[Abstract/Free Full Text]
• OPTICAN, L. M., RICHMOND, B. J. Temporal encoding of two-dimensional patterns by single units in primate inferior temporal cortex. III. Information theoretic analysis. J. Neurophysiol. 57: 162-178, 1987.[Abstract/Free Full Text]
• PARADISO, M. A. A theory for the use of visual orientation information which exploits the columnar structure of striate cortex. Biol. Cybern. 58: 35-49, 1988.[Medline]
• PELLIONISZ, A., LLINÁS, R. Tensor network theory of the metaorganization of functional geometrics in the central nervous system. Neuroscience 16: 245-273, 1985.[Medline]
• PELLIONISZ, A. J. Coordination: a vector-matrix description of transformation of overcomplete CNS coordinates and a tensorial solution using the Moore-Penrose generalized inverse. J. Theor. Biol. 110: 353-375, 1984.[Medline]
• POGGIO, T. A theory of how the brain might work. In: Cold Spring Harbor Symposia Quantitative Biology, The Brain., . Plainview, NY: Cold Spring Harbor Laboratory Press, 1990, vol. LV, p. 899-910
• POGGIO, T., GIROSI, F. Regularization algorithms for learning that are equivalent to multilayer networks. Science 247: 978-982, 1990.[Abstract/Free Full Text]
• POGGIO, T., TORRE, V., KOCH, C. Computational vision and regularization theory. Nature 317: 314-319, 1985.[Medline]
• POUGET, A., SEJNOWSKI, T. J. A neural model of the cortical representation of egocentric distance. Cereb. Cortex 4: 314-329, 1994.
• POUGET, A., SEJNOWSKI, T. J. Spatial transformations in the parietal cortex using basis functions. J. Cog. Neurosci. 9: 222-237, 1997.
• POUGET, A., ZHANG, K.-C. Statistically efficient estimations using cortical lateral connections. In: Advances in Neural Information Processing Systems, edited by M. C. Mozer, M. I. Jordan, and T. Petsche . Cambridge, MA: MIT Press, 1997, vol. 9, p. 97-103
• RANCK, J. B. JR Studies on single neurons in dorsal hippocampal formation and septum in unrestrained rats. I. Behavioral correlates and firing repertoires. Exp. Neurol. 41: 461-555, 1973.[Medline]
• RAO, C. R. In: Linear Statistical Inference and Its Applications., . New York: Wiley, 1965
• REDISH, A. D., TOURETZKY, D. S. Modeling interactions of the rat's place and head direction system. In: Advances in Neural Information Processing Systems, edited by D. S. Touretzky, M. C. Mozer, and M. E. Hasselmo . Cambridge, MA: MIT Press, 1996, vol. 8, p. 61-67
• RIEKE, F., WARLAND, D., DE RUYTER VAN STEVENINCK, R., BIALEK, W. In: Spikes: Exploring the Neural Code., . Cambridge, MA: MIT Press, 1997
• ROLLS, E. T., ROBERTSON, R. G., GEORGES-FRANCOIS, P. The representation of space in the primate hippocampus. Soc. Neurosci. Abstr. 21: 1494, 1995.
• ROTENBERG, A., MAYFORD, M., HAWKINS, R. D., KANDEL, E. R., MULLER, R. U. Mice expressing activated CaMKII lack low frequency LTP and do not form stable place cells in the CA1 region of the hippocampus. Cell 87: 1351-1361, 1996.[Medline]
• SALINAS, E., ABBOTT, L. F. Vector reconstruction from firing rates. J. Comput. Neurosci. 1: 89-107, 1994.[Medline]
• SAMSONOVICH, A., MCNAUGHTON, B. L. Attractor-map-based path integrator model of the hippocampus reproduces the phase precession phenomenon. Soc. Neurosci. Abstr. 22: 1872, 1996.
• SAMSONOVICH, A., MCNAUGHTON, B. L. Path integration and cognitive mapping in a continuous attractor neural network model. J. Neurosci. 17: 5900-5920, 1997.[Abstract/Free Full Text]
• SANGER, T. D. Theoretical considerations for the analysis of population coding in motor cortex. Neural Comput. 6: 12-21, 1994.
• SANGER, T. D. Probability density estimation for the interpretation of neural population codes. J. Neurophysiol. 76: 2790-2793, 1996.[Abstract/Free Full Text]
• SCHARF, L. L. In: Statistical Signal Processing: Detection, Estimation, and Time Series Analysis., . Reading, MA: Addison-Wesley, 1991
• SCHWARTZ, A. B. Motor cortical activity during drawing movements: Population representation during sinusoid tracing. J. Neurophysiol. 70: 28-36, 1993.[Abstract/Free Full Text]
• SCHWARTZ, A. B. Direct cortical representation of drawing. Science 265: 540-542, 1994.[Abstract/Free Full Text]
• SEUNG, H. S., SOMPOLINSKY, H. Simple models for reading neuronal population codes. Proc. Natl. Acad. Sci. USA 90: 10749-10753, 1993.[Abstract/Free Full Text]
• SKAGGS, W. E., MCNAUGHTON, B. L. Replay of neuronal firing sequences in rat hippocampus during sleep following spatial experience. Science 271: 1870-1873, 1996.[Abstract]
• SKAGGS, W. E., MCNAUGHTON, B. L., GOTHARD, K. M., MARKUS, E. J. An information-theoretic approach to deciphering the hippocampal code. In: Advances in Neural Information Processing Systems, edited by S. J. Hanson, J. D. Cowan, and C. L. Giles . Cambridge, MA: MIT Press, 1993, vol. 5, p. 1030-1037
• SKAGGS, W. E., MCNAUGHTON, B. L., WILSON, M. A., BARNES, C. A. Theta phase precession in hippocampal neuronal populations and the compression of temporal sequences. Hippocampus 6: 149-172, 1996.[Medline]
• SNIPPE, H. P. Parameter extraction from population codes: a critical assessment. Neural Comput. 8: 511-529.
• SOMERS, D. C., NELSON, S. B., SUR, M. An emergent model of orientation selectivity in cat visual cortical simple cells. J. Neurosci. 15: 5448-5465, 1995.[Abstract]
• STEWART, M., FOX, S. E. Do septal neurons pace the hippocampal theta rhythm? Trends Neurosci. 13: 163-168, 1990.[Medline]
• TAUBE, J. S., MULLER, R. U. Head direction cell activity in the anterior thalamic nuclei, but not the postsubiculum, predicts the animal's future directional heading. Soc. Neurosci. Abstr. 21: 946, 1995.
• TAUBE, J. S., MULLER, R. U., RANCK, J. B. JR Head-direction cells recorded from the postsubiculum in freely moving rats. I. Description and quantitative analysis. J. Neurosci. 10: 420-435, 1990.[Abstract]
• TIKHONOV, A. N., ARSENIN, V. Y. In: Solutions of Ill-Posed Problems., . New York: Wiley, 1977
• TREVES, A., BARNES, C. A., ROLLS, E. T. Quantitative analysis of network models and of hippocampus data. In: Perception, Memory, and Emotion: Frontiers in Neuroscience, edited by T. Ono, B. L. McNaughton, S. Molotchinkoff, E. T. Rolls, and H. Nishijo . Oxford, UK: Pergamon, 1996, p. 567-579
• TSODYKS, M. V., SEJNOWSKI, T. J. Associative memory and hippocampal place cells. Int. J. Neural Syst. 6, Suppl.: 81-86, 1995a.
• TSODYKS, M. V., SEJNOWSKI, T. J. Rapid state switching in balanced cortical network model. Network 6: 111-124, 1995b.
• TSODYKS, M. V., SKAGGS, W. E., SEJNOWSKI, T. J., MCNAUGHTON, B. L. Population dynamics and theta rhythm phase precession of hippocampal place cell firing: a spiking neuron model. Hippocampus 6: 271-180, 1996.[Medline]
• TUCKWELL, H. C. In: Introduction to Theoretical Neurobiology: Nonlinear and Stochastic Theories., . Cambridge, UK: Cambridge Univ. Press, 1988, vol. 2
• ULLMAN, S., BASRI, R. Recognition by linear combinations of models. IEEE Trans. Pat. Anal. Mach. Intel. 13: 992-1006, 1991.
• VANDERWOLF, C. H., KRAMIS, R., GILLESPIE, L. A., BLAND, B. H. Hippocampal rhythmical slow activity and neocortical low voltage fast activity: relations to behavior. In: The Hippocampus: Neurophysiology and Behavior, edited by R. L. Isaacson, and K. H. Pribram . New York: Plenum, 1975, vol. 2, p. 101-128
• (corrections in vol. 264, p. 16) WILSON, M. A., MCNAUGHTON, B. L. Dynamics of the hippocampal ensemble code for space. Science 261: 1055-1058, 1993.[Abstract/Free Full Text]
• WILSON, M. A., MCNAUGHTON, B. L. Reactivation of hippocampal ensemble memories during sleep. Science 265: 676-679, 1994.[Abstract/Free Full Text]
• ZACKS, S. In: The Theory of Statistical Inference., . New York: Wiley, 1971
• ZEMEL, R. S., DAYAN, P., POUGET, A. Probabilistic interpretation of population codes. In: Advances in Neural Information Processing Systems, edited by M. C. Mozer, M. I. Jordan, and T. Petsche . Cambridge, MA: MIT Press, 1997, vol. 9, p. 676-682
• ZHANG, K.-C. Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. J. Neurosci. 16: 2112-2126, 1996.[Abstract/Free Full Text]
• ZIPSER, D. A computational model of hippocampal place fields. Behav. Neurosci. 99: 1006-1018, 1985.[Medline]
• ZIPSER, D., ANDERSEN, R. A. A back-propagation programmed network that simulates response properties of a subset of posterior parietal neurons. Nature 331: 679-684, 1988.[Medline]
• ZOHARY, E. Population coding of visual stimuli by cortical neurons tuned to more than one dimension. Biol. Cybern. 66: 265-272, 1992.[Medline]
• ZOHARY, E., SHADLEN, M. N., NEWSOME, W. T. Correlated neuronal discharge rate and its implications for psychophysical performance. Nature 370: 140-143, 1994.[Medline]

0022-3077/98 $5.00 Copyright ©1998 The American Physiological Society