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