Abstract
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 multielectrode 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 informationtheoretic 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
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érRao 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).
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 singletrial reconstruction, phenomena such as erratic jumps in reconstructed trajectory easily can become apparent (Fig. 9).
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 figure8 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 figure8 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 3min 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
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
The firing rate distribution f(x) can be interpreted in terms of conditional probability as follows
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
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.
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.
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
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 MoorePenrose 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
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
To find the optimal template functions g_{i}
(x), we use vectormatrix notation to rewrite Eq. 10
as
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 = UD
^{†}
V
^{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
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
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
In summary, given the numbers of spikes (n
_{1}, n
_{2}, . . . , n
_{N}) fired by all the cells during a given time interval, combine the chosen basis functions φ_{i}(x) linearly as
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
_{1} or f
_{2}, 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
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
_{1} and f
_{2} as column vectors and define F = [f
_{1}, f
_{2}] as a 2 × 2 matrix. Equation 22
is equivalent to
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
_{1}(x), f
_{2}(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
_{1}, n
_{2}, . . . , 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 =
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
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 x̂
_{vector} in Eq. 28
is identical to the peak position of the following scalar function
POPULATION VECTOR WITH SCALING.
For place cells, the population vector needs to be scaled by the total activity so that the reconstructed position
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 2D 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
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 informationtheoretical 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 stimulusresponse relation can be obtained empirically as a lookup table (Földiák 1993).
Basic method
Assume that using the methods in the preceding section we have measured the firing rate maps f _{1}(x), f _{2}(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
_{1}, n
_{2}, . . . , 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
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
The final formula obtained from Eq. 34
reads
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
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
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
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 onestep 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
_{t−1} 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
_{t−1}) from the general Bayes formula P(x
_{t−1}‖x
_{t}, n
_{t})P(x
_{t}, n
_{t}) = P(x
_{t}‖n
_{t}, x
_{t−1})P(n
_{t}, x
_{t−1}) and the assumption that
Now the reconstructed position at time t is
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
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
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}/2σ^{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
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 2D 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.^{1}
TIME WINDOW.
A time window with a duration of ∼1 s was a good choice for the Bayesian twostep method for these data (Fig. 6 A). 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 twostep 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 6 B 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 selfreconstruction because all parameters for the reconstruction algorithm were determined using data from only the first half session.
Upon closer examination of Fig. 6 B, 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 selfreconstruction curve (light symbols). To quantify the drift, each curve was fitted with a fifthorder 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. 6 B, 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 selfreconstruction (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.
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 3min 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 3min 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 figure8 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 twostep 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 ). 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érRao 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 2D 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 2D Gaussian
Here we only state the final results of the detailed derivation provided in . The simplest result is that reconstruction error cannot be smaller than
Another equivalent error formula is
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 ().
Unbiased reconstruction is assumed implicitly when the CramérRao 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
The minimal reconstruction error can be reached by the Bayesian method (cf. ). 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
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 1D situation that corresponds to a 2D 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
In the same situation, the average angular error of the population vector method also can be estimated directly (). The result is
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 twostep 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 . For animal 1, the average width of place fields was 〈σ〉 = 10.1 cm and
Because the Bayesian twostep 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 onestep 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érRao lower bound (Figs. 10 and 12 and ). 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érRao 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 (). 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 twostep method, the empirical value was 0.8705 for animal 1 and 0.8855 for animal 2, surprisingly close to the value of the 2D correction factor ℱ_{2} = 0.8862 but quite different from ℱ_{1} and ℱ_{3} (Table B1). However, the actual error distribution differed from the theoretical probability density ρ_{2}(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 2D 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 longterm 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
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 1D space. Examples of 1D tuning curves include the orientation tuning in visual cortex and directional tuning of headdirection cells (Taube et al. 1990). In contrast, in3D space or space of higher dimensions, a sharper tuning width implies worse information encoding. So the 2D 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 x̂ in all reconstruction methods considered above can be expressed as
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
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
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
Indeed, linear combination of basis functions is biologically plausible computational strategy that has been used by many authors, such as in 3D 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 winnertakeall operation on the activity distribution in the second layer after the first feedforward step. An exact winnertakeall 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 onestep Bayesian method can be implemented simply as constant inputs to the second layer. To implement the additional bias for the twostep Bayesian method, a type of shortterm memory is needed for the winnertakeall 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 winnertakeall, 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 (BenYishai 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érRao lowerbound 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érRao 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 for1D 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 2D space, the variable could as well be discrete and disjoint, which might be more suitable for representing distinct classes of objects or categories.
Acknowledgments
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 NS20331.
Footnotes

Address reprint requests to K. Zhang.

1 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 loglog 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.
 Copyright © 1998 the American Physiological Society
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
The reciprocal basis functions obtained by Eq. EA2
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
These formulas are related. In Eq. EA5 , the D term is related to firing rate variability (see below). For large time window τ, the D term vanishes, and Eq. EA5 is reduced to Eq. EA4 . When the visiting probability P_{k} is the same for all positions, Eq. EA4 is reduced to the original reciprocal basis Eq. EA2 .
Equation EA5
has been derived by minimizing the error function
When firing rate variability is ignored, the equation
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érRao inequality.
As an intuitive example, consider a 1D 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
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):
The squared error for any unbiased estimator of position x cannot exceed the CramérRao bound, namely
When more than one cell is used in reconstruction, as long as they are independent, the total Fisher information is the sum
So far, we have only considered a 1D problem. Extension to 2D 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
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 Ddimensional space
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
_{1} for all N cells should be
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
The second method uses the variance
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. EB17
, which gives
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 2D 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 2D distribution f(i, j) to obtain a 1D distribution f(i) and then assume a 1D Gaussian distribution to obtain the 1D tuning width
Distribution of reconstruction error and correction factor
A correction factor is needed because the CramérRao bound only gives an estimate of the variance of unbiased reconstruction, that is, the average square of the reconstruction error 〈r
^{2}〉. Its square root
For reconstruction in Ddimensional space, it is natural to assume that the distribution of the error vector obeys a Ddimensional Gaussian distribution
The probability density ρ_{D}(r) for reconstruction error can thus be obtained from Eqs. B22
and B23
:
The correction factor, defined as 〈r〉/
The results for 1, 2, and 3D 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 Ddimensional 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).
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érRao 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érRao bound in the limit of large sample (Cramér 1946).
Consider the Fisher information for a 1D reconstruction under a uniform prior probability:
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
APPENDIX D: ERROR OF POPULATION VECTOR UNDER COSINE TUNING
2D 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 xaxis, u = (1, 0). The average firing rate of cell i is
Because the test direction is assumed to be along the xaxis, population vector U is also close to the xaxis. Thus ‖U‖ ≈ U_{x}
and the angular error of the population vector can be approximately as
3D work space
We need to consider only a single test direction along the xaxis: u = (1, 0, 0). The method is very similar to that in the case of 2D work space, except that now