Abstract
The firing of hippocampal principal cells in freely running rats exhibits a progressive phase retardation as the animal passes through a cell's “place” field. This “phase precession” is more complex than a simple linear shift of phase with position. In the present paper, phase precession is quantitatively analyzed by fitting multiple (1–3) normal probability density functions to the phase versus position distribution of spikes in rats making repeated traversals of the place fields. The parameters were estimated by the Expectation Maximization method. Three data sets including CA1 and DG place cells were analyzed. Although the phaseposition distributions vary among different cells and regions, this complexity is well described by a superposition of two normal distribution functions, suggesting that the firing behavior consists of two components. This conclusion is supported by the existence of two distinct maxima in the mean spike density in the phase versus position plane. In one component, firing phase shifts over a range of about 180°. The second component, which occurs near the end of the traversal of the place field, exhibits a low correlation between phase and position and is antiphase with the phaseshift component. The functional implications of the two components are discussed with respect to their possible contribution to learning and memory mechanisms.
INTRODUCTION
Since the 1970s' the hippocampal cognitive map theory has been developed on the basis of the assumption that the firing rate of individual cells encodes the animal's position (O'Keefe and Dostrovsky 1971;O'Keefe and Nadel 1978). O'Keefe and Recce (1993), however, showed that spatially specific firing of pyramidal cells in the CA1 region of the rat hippocampus also has a characteristic temporal relationship with the 7Hz theta field potential recorded near the hippocampal fissure. As the rat enters the place field of a given cell, the initial spikes occur late in the theta cycle. As the rat passes through the place field, the phase advances progressively. The last spikes fired by the cell as the rat leaves the place field occur almost 360° earlier in the theta cycle, and, for unimodal place fields, phase never shifts more than 360° over the range of the field.
O'Keefe and Recce (1993) called this phenomenon “theta phase precession.” Using data from animals shuttling back and forth on a linear track, and considering the two directions of travel independently, they applied linear regression analyses to the theta phase versus position plot of individual units and concluded that firing phase is highly correlated with spatial location. Thus within a given trajectory on a linear track, firing phase provides substantial additional information about the animal's location beyond that provided by firing rate.
Skaggs et al. (1996) examined the phase precession phenomenon using data from largescale neuronal ensemble recordings in rats running unidirectionally around a triangular track. They pointed out that two consequences of the phase shift phenomenon were that repeating, overlapping segments of the global sequence of place fields are expressed in compressed form within each theta cycle, and that as a result, the neural population codes at the beginnings and ends of a theta cycle are essentially uncorrelated. They proposed that this compression and decorrelation could enhance the learning of sequences by means of Hebbian longterm potentiation (LTP), which has a time constant of <50 ms (Bi and Poo 1998; Gustaffson and Wigström 1986; Larson and Lynch 1989; Levy and Steward 1983).
Therefore it is important to clarify the neural mechanism generating theta phase precession to understand how this process might facilitate learning and memory of behavioral sequences. In this regard, a number of theoretical models have been proposed (Bose et al. 2000; Jensen and Lisman 1996b; O'Keefe and Recce 1993; Tsodyks et al. 1996;Wallenstein and Hasselmo 1997). All of these models are based on the assumption that the phase shift is smooth and ranges up to one period of the theta cycle, as reported by O'Keefe and Recce (1993). On the contrary, Skaggs et al. (1996)observed that the relation between phase and location is typically not a simple linear shift, but has more complex dynamics, which may vary from cell to cell or from place field to place field for a given cell. They showed several types of complexity. In some cases the phase advance accelerates, resulting in a downward curve of the phase versus position plot, and in other cases the shift was not smooth, but exhibited two or more distinct clusters in the theta phase versus position plots. Examples are shown in Fig.1. From this viewpoint, the example shown in Fig. 3 of O'Keefe and Recce (1993) is a clear case of two apparently distinct clusters. Further evidence of such complexity is seen in Fig. 3B of Skaggs et al. (1996), in which the average population activity distribution versus theta phase profile exhibits bimodality, and in their Fig. 6,B–D, which shows bimodal phase distributions in the later portions of individual place fields.
Such observations suggest that phase precession is a more complex physiological phenomenon than previously hypothesized. To understand the complexity of the firing phase behavior, it is important to develop some form of quantitative description that captures the possible bimodality of the firingphase distributions. In the present paper, we hypothesize that there are two components in the form of normal probability density functions, in which the firing phase behavior consists of two different kinds of neural activity, and in which further complexity of the firing phase behavior may possibly appear as a consequence of variable contributions of the two components in different cases. The following sections are therefore devoted to the quantitative description of the firing phase behavior based on the twocomponent hypothesis. The results are submitted to further discussion on the possible neural mechanism of phase precession and its functional role in learning and memory.
METHODS
Electrophysiological data obtained during three different behavioral experiments were used in the present investigation.Data set 1 is identical to that reported in Skaggs et al. (1996). In that experiment the rats were required to run repeatedly around a triangular track maze (80 cm per side) for food reinforcement delivered at the midpoint of each side. Data set 2 is from a similar experiment in which the apparatus consisted of a circular track maze 120 cm diame. Data set 3 is from animals shuttling back and forth on a simple linear track 180 cm in length. The motion of the rat was always counterclockwise in data sets 1 and 2. In data sets 1 and3, cells were recorded from simultaneously in CA1 and the dentate gyrus (DG). For position analysis, lightemitting diodes mounted on the rats' headstages were tracked at 60 Hz by a video system that was synchronized with the electrophysiological signal time stamps. Spike trains from multiple hippocampal pyramidal or granule cells were recorded using multiple tetrode probes as described bySkaggs et al. (1996). Putative single neurons were isolated using an interactive clustering procedure on the spike waveform and amplitude parameters. After spike sorting, suitable spike trains were selected for phase precession analysis according to the following criteria: 1) a sufficient number of spikes (a minimum of 200 was necessary to ensure convergence in the statistical analyses described below); 2) an approximately unimodal distribution of spikes on the position coordinate (place field); and3) the place field being not too close to the feeding position. Criterion 3 was applied for two reasons. First, it essentially eliminated spikes that occurred while the animal was not in a state of forward locomotion. Second, because the analysis involved collapsing the position coordinate to a single dimension, it was necessary to eliminate spikes that occurred while the rat was turning around or moving its head from side to side, which frequently occurred near the reinforcement sites. The analyses below are applied to a number of spike trains, some of which expressed more than one place field. For simplicity, each place field is considered separately and the term “unit” is used here to refer to a place field. Place fields were defined subjectively as approximately unimodal regions of elevated firing spanning a distance of about 30 cm.
Data set 1 consists of 29 units from 26 CA1 cells and 4 units from 4 DG cells from a single recording session with 145 laps.Data set 2 consists of 15 units from 12 CA1 cells from a single session with 14 laps. Data set 3 consists of 10 units from 7 CA1 cells and 20 units from 13 DG cells from 4 sessions with 12–65 laps. The differences in the number of laps among data sets contributes to differences in the total number of spikes per unit. Thus the analysis focuses on data set 1, which contained the most laps, whereas, results from data sets 2 and 3 are likely to be less quantitatively precise and are included primarily for qualitative confirmation of the essential findings.
For theta phase analysis, electroencephalographs (EEGs) from an electrode at or near the hippocampal fissure were filtered in the range of 6–10 Hz (see details of experimental method in Skaggs et al. 1996). The time stamp of the positive peak of the filtered wave was used as the zero phase reference for the theta rhythm. Note, however, that differences in electrode location contribute to apparent phase shifts in the theta rhythm, and so the zero point should be taken as a relative reference in each session.
The theta phase of the ith spike of the sth unit is defined by
The rat position p
Estimation of the probability density function by means of the expectation maximization (EM) algorithm
For the purpose of analysis, the spike distribution in thepφ plane was assumed to conform to a set of normal probability density functions. The EM algorithm (Dempster et al. 1977; Render and Walker 1984) was applied to estimate the parameters of these functions. The EM algorithm is sometimes used for cluster analysis and is also identical with maximum likelihood estimation of incomplete data. This algorithm is defined briefly as follows. The probability density of a multivariate space is assumed to be described by the superposition of a number of normal density functions. The parameters of the individual normal density functions are initially selected at random and are changed iteratively so that the expectation of the observed spike distribution increases to a maximum. When the values of the parameters converge, the result is the bestfitting composite probability density function for the given number of components.
In the simplest case of a single normal density function, the probability density function in the (x, y) plane is represented by
When the probability density function is assumed to consist ofN normal density functions, the functionP(x, y) becomes
Before the EM algorithm is applied to the spike data of place units, a complicating factor in the present application must be considered. Since the variable φ is a modulus, the twodimensional plane defined by the pφ coordinate is actually a cylinder, whereas the normal density function is defined in a plane and does not satisfy the periodic boundary condition. A related difficulty appears when one cuts the cylinder at any arbitrary φ value. Even if a single cluster structure is obvious in the cylindrical coordinate, the distribution in the planar representation would generally consist of two clusters because of the artificial boundary effect, unless the cut happened to be exactly at the minimum of the distribution. For two or more clusters with variable peak locations in the φ coordinate, spurious splitting of clusters is unavoidable.
To minimize the above difficulties, the EM algorithm was applied to distributions in which the φ coordinate was replicated over multiple theta periods. When the number of cycles is M, andL normal density functions are assumed in a cycle, the EM algorithm is applied with n = LM inEq. 3a . When M is a sufficiently large number, the parameters of the components in the middle parts tend to give essentially periodic distributions, free from the boundary effects (see example in Fig. 2.).
The new variable Φ representing M periods of theta phase is introduced below. Φ is not a modulus but a real number defined in the interval [φ_{0}, φ_{0} +M]. The value φ_{0} is usually chosen as the theta phase value at the trough of the population firing distribution in each data set, because this defines the natural boundary with the minimal number of clusters in the pφ plane (i.e., the boundary that does not cause obvious artificial clusters). In the following, M = 5 periods of theta cycles is used. The value of L is varied from 1 to 3 for comparison. Furthermore, the results are evaluated with the following criteria for periodicity.
 1)
 The arrangement of normal density functions conform to periodic density functions, except for the first or last cycles in the neighborhood of the boundary.
 2)
 The estimated individual normal functions are effectively localized within each theta period (recall that the phase shift for a place field spans only a single cycle).
Conditions 1) and 2) were first evaluated by eye. If the two conditions were obviously violated, the unit was removed from the statistics below. Condition 1) was considered first at the start of the EM algorithm as follows. The initial values of parameters in Eq. 3a were chosen in a manner consistent with the periodicity along the Φ axis: the initial location parameters (means) in Φ and P are random values within each theta cycle. The initial values of variance in the Φ axis are 0.01. The initial values of variance in the P axis are 0.6 (data set 1), 0.43 (data set 2) and 0.78 (data set 3). The initial values of covariance are 0, and the initial values of α_{i} are 1/N. After the maximization algorithm was applied, the periodicity of the resultant distribution was checked in the range from the second to the fourth period, i.e., [φ_{0} + 1, φ_{0} + 3] in Φ. If the periodicity constraint was obviously violated, the initial values of the means were changed randomly for the next trial of the EM algorithm. In the present case, such intervention was unnecessary except in a few cases. Most units yielded periodicity after the first trial. Conditions 1 and2 were satisfied in each data set containing a minimum of 200 spikes. In the following sections, the theta phase is described by either φ or Φ. Note that φ is modulus 1 and Φ is a real variable in the range [0, 5].
To illustrate a typical case of the above procedure, Fig. 2,A and B, shows an example of the spike distribution in the pφ and pΦ planes. The distribution of spikes is sparse around φ_{0}
= 0.25. The result of applying the EM algorithm, withL = 2, to this unit is shown in Fig. 2, Cand D. Figure 2
C shows the 10 individual normal density functions, each of which is denoted by a contour line at the level of 1/e. Figure 2
D shows the contour plot of the estimated probability density function composed of 10 normal density functions. It is seen that each normal density function is localized essentially within a single theta period. Individual functions are roughly periodically arranged along the Φ axis. The clusters near the upper and lower boundaries are different in detail from the others nearer the middle. Thus the result in the middle cycle is essentially free from the boundary effect and represents a reasonable approximation of the probability density function in a cylindrical coordinate system. The L normal density functions whose centers are located at the third theta cycle were taken as the final estimated parameters in the pφ cylinder coordinate. The relative weights of individual functions {α*_{j}} are defined by normalization among {α_{j}
} forL functions. Thus the estimated parameters with a given value of L are obtained for each unit as follows
The validity of the parameter estimation can be evaluated by the log likelihood function (McLachlan and Whiten 1996;Render and Walker 1984). The parameter λ_{L} is defined by the log likelihood function with a conditional probability as follows
RESULTS
General results
Several examples of results obtained in data set 1 by the EM algorithm are shown in Fig.3. The results withL = 1, 2, and 3 are shown for the sake of comparison. In Fig. 3 A, the result of L = 1 gives a probability density function with a single peak. The distribution with the shape of an oblique ellipse could correspond to the linear phase shift in one theta cycle. The principal axis of the ellipse (i.e., the longer one of the 2 orthogonal axes of the ellipse) corresponds to the approximated linear regression line. When L = 2 is used, two normal density functions, A and B, are found to have almost equal weights. The resultant probability density function shows two peaks reflecting the peaks of individual normal density functions. The normal density function termed A in the figure has rather smaller variance and a higher peak, than B. When the ridge of the distribution over one theta cycle is traced across the two peaks, an “Sshaped” form is apparent. It also appears as two distinct clusters. In the case of L = 3, the resultant probability density function is found to show two peaks. The two peaks with L = 3 appear similar to the two peaks with L = 2. That is, the third normal density function is apparently not significant but merely describes random details of the sample.
Figure 3 B shows another example. There is an apparent acceleration of the phase shift along the position coordinate of the place field with two peaks in the probability density function with either L = 2 or L = 3. Thus two examples in Fig. 3, A and B, have some differences in their apparent shape but show probability density functions with two peaks with both L = 2 andL =3.
Figure 3 C shows yet another type of unit. The probability density function shows a single peak for L = 1, 2, and 3. Thus the shape of the single peak appears robustly in this unit. The shape and location of the distribution is similar to the upper cluster in the case with L = 2 or 3 in Fig. 3, A orB.
Table 1 shows a summary of the analyses of all units of data sets 1–3. Convergence of the parameter estimation by the EM algorithm was established in all units shown in the table. It is seen that the probability density functions forL = 1 ∼ 3 have single or double peak(s) in 26/29 units of data set 1, and in 9/15 units of data set 2. Thus in these units, the qualitative shape of the probability density function obtained for L = 1 or 2 is stable for further increases in L. In data set 3, the probability density function with a single peak is kept forL = 2 in a half of units. The result forL = 3 of data set 3 is not shown because the small number of spikes causes difficulty in the convergence of the parameters.
The parameters λ_{L} with various values for individual units of data sets 1 and 2 are shown in Fig. 4. It is observed that the difference between λ_{1} and λ_{2} is larger than that between λ_{2} and λ_{3} in every unit. This means that the proportion of variance in the raw data that is explained by the fitted probability density functions saturates in the range of L < 3. It should be noted that even in the case of a single cluster with a single peak for all L,the value of λ_{L} will increase monotonically because extraneous normal functions are fitted to the details of the spike distribution. Therefore an increase in λ_{L} cannot be used by itself to infer the existence of multiple clusters in the data; however, the saturation tendency indicates that further analyses beyond L = 3 are unlikely to be meaningful. In the following sections, we focus on the results for L = 2, for quantitative evaluation of spike distribution.
Probability density functions of CA1 units
In the case of L = 2, we distinguish the two normal density functions of each unit by the numbers N1 and N2 or by suffixes 1 and 2. The number 1 is used for the normal density function whose peak value on the p axis is smaller (i.e., near the entry of the place field), so that μ_{p1} < μ_{p2}. When the relative weights of the two density functions are compared over all units, the mean of α
Figure 6 summarizes the results of CA1 units in data sets 1and 2. The relative locations of the two normal density functions are summarized in terms of the center of the distribution in Fig. 5, A and B. The center of the normal function N1 is used as the origin of the abscissa by using (μ_{pi} − μ_{p1},μ _{φi}) (i = 1, 2). The average of the distance in between N1 and N2 along the position axis, μ_{p2} −μ _{p1} is 8.86 ± 3.27 cm indata set 1 and 12.98 ± 4.30 cm in data set 2. In data set 1 shown in Fig. 5 A, the distributions of μ_{φ1} and μ_{φ2} on the ordinate are almost separated from each other. The center of mass of N1 on the ordinate is 2.94 ± 0.08 in Φ, i.e., about 0 in φ, and the center of N2 is around 2.55 ± 0.15 in Φ, i.e., 0.55. That is, the two are almost antiphase to each other. The antiphase relation between N1 and N2 is similarly observed in data set 2 as shown in Fig.5 B.
The density of N1 and N2 in individual units is shown in Fig.6. In the spatial domain, the activity of every unit is initiated by N1 and terminated by N2, while the two normal distributions are significantly overlapped with each other. It is seen that N1 has a dominant density around the upper half theta cycle. On the other hand, N2 mostly distributes at the lower half. Thus N1 and N2 are characterized by their relative distributions within the place field and in theta phase. That is, the two normal density functions are typically identified as one within the initial portion of the field and the other toward the end of the field, with almost antiphase relations to each other. Since the phase relation is relative, the location within the place field is conveniently used for the classification of N1 and N2.
The relation between position and theta phase in N1 and N2 are summarized in the lower part of Fig. 5. The correlation coefficientsr _{pφj} in N1 and N2 are shown in Fig. 5, C and D. It is obvious that N1 has high negative correlation of phase and location in both of data sets 1 and 2. The negative value means that, within N1, the phase of unit firing advances as the rat traverses the place field. On the other hand, the correlation coefficient of N2 has zero mean and large variance. Thus N1 essentially represents what is called theta phase precession or the phase shift component, whereas N2 appears to endow complexity and variety in firing phase behavior but does not entail a systematic phase shift. The gradient of the principal axis of N1 is plotted in Fig. 5, E and F, which shows that the gradient distributes around −0.04 ± 0.03 theta cycle/cm in data set 1 and −0.02 ± 0.01 in data set 2. As is shown in the next section, the range of the distribution of CA1 units is also similar to that of DG units in the same data set.
Figure 7 shows the probability density functions obtained as the average of the parameters {Π_{2}} for individual units in data set 1. Note that position of the center of N1 or N2 is averaged by using relative location, (μ_{pi} − μ_{p1}, μ_{φi}) (i = 1, 2). It is seen that the average of all CA1 units in A shows the single peak. The components derived from N1 and N2 are mixed to a single cluster. At a given value of relative position, the probability density function has a single maximum value of φ. The maximum point shifts downward with relative location, in a gradually accelerating fashion. The initial linear phase shift comes from N1 and the late acceleration comes from N2. As a whole, the distribution covers almost a full theta cycle.
When the average is performed over units with double peaks (n = 23) as shown in B, a second peak appears late in the place field, antiphase to the first peak. The average distribution of units with a single peak (n = 6) is shown in C, where N1 and N2 appear to represent a single cluster in the upper half cycle very similar to the peak in the initial half of the place field in B. Thus grouping units according to the number of peaks provides a qualitative classification scheme. The two groups differ in the presence or absence of the additive peak caused by an increased or decreased N2 contribution. It should be noted that the difference in the number of peaks does not necessarily imply distinct properties of N1 and N2. A gradual change of any parameter can result in a sudden change in the number of peaks as a bifurcation. Further analysis shows that the probability density function in Fig. 7 A is close to the bifurcation point and sensitive to small parameter change, and therefore the subtraction of six units with a single peak from the average ensemble caused a change from a single peak into a double one. Nevertheless, the number of peaks provides a means of classifying units with respect to the contribution of N2.
Similar properties in the relation between the two peaks are also found in data set 2 as shown in Fig.8.
Probability density functions of DG units
For DG units, two normal density functions N1 and N2 were defined according to the sequence of their positions, as in the previous section. The estimated parameters are summarized in Fig.9. For the sake of comparison, the result of CA1 units in data set 3 is also shown in Fig.10. As is shown in Fig. 9, Aand B, the relative location of the two peaks appears similar to that of CA1 units. The range of the distribution of N2 along Φ axis is, however, a little smaller than that in CA1. As shown in Fig. 9, C and D, the correlation coefficientr of N1 in DG units is strongly negative as in CA1 units. The gradient of the principal axis of N1 of DG units is shown in Fig.9, E and F. The ranges of the distributions of the gradients in CA1 and DG are similar in each data set.
Figure 11 shows the probability density functions obtained as an average of the parameters {Π_{2}}, for DG units. The average indata set 1 (Fig. 11 A) shows a clear single peak naturally because all four DG units in this set have a single peak. It should be noted that the shape of the distribution is very similar to that obtained for CA1 units with a single peak in data set 1. The difference between the two is mostly restricted to the phase shift.
Since onehalf of DG units in data set 3 show double peaks as shown in Table 1, the parameters of the probability density function are averaged over all units in Fig. 11 B, over units with double peaks in Fig. 11 C and over units with a single peak in Fig. 11 D. It is seen that all these results show very similar distributions with a single peak. The second peak of individual units disappeared in Fig. 11 C after averaging, because the N2 components are randomly located around the peak of N1. That is, the collective distribution of N2 does not result in a second peak in the average. Rather the profiles remain single peaked as in data set 1.
Comparison of collective firing properties obtained in various ways
In the preceding, the collective distributions for some units show double peaks; however, it is possible that this effect might be an artifact of the description by two normal distributions. In this section, the collective firing properties of individual units are further evaluated in other ways to avoid this potential artifact, and provide further evidence for the existence of double peaks.
First, the probability density functions of individual units were summed over a given set of units in two ways: after first aligning the data on the center of N1 or after aligning the data on the center of N2. Putting the number of units n and L = 2 in the EM algorithm, the resultant function is given by summation of 2n normal density functions. That is, the resultant function generally may become more complex and the number of peaks may be either less or more than two. Figure 12 shows the averages for CA1 and DG units in data set 1 with the center of N1 as the origin. It is seen that results for CA1 and DG units are essentially the same as those of the probability density function with averaged parameters in Figs. 7 and 11 A. They have in common an initial region of phase shift over the first half of the theta cycle. The net phase shift in the early component is similar in DG and CA1 (about 0.2). The result obtained from CA1 units with double peaks (Fig. 12 B) shows the additional peak at the second half cycle. Figure 13 shows the results obtained by summation with data aligned on the center of N2. The resultant function tends to emphasize the peak derived from N2, but the difference between Figs. 12 and 13 is small. The general form of the average functions thus does not depend much on whether the data are aligned on N1 or N2. Both are quite similar to the results obtained by averaging the individual parameters of the probability density functions, indicating that the collective properties shown above are relatively consistent over the population of units.
For the sake of comparison, Fig. 14shows the actual average density distribution of spikes of CA1 and DG units in data set 1. It is observed that the actual average density plots exhibit roughly the same patterns of contour lines as the EM analysis. That is, CA1 units show single or double peaks and DG units consistently show a single peak as observed in the EM analysis with L = 2. It should be noted that number of peaks in the density plot depends inversely on the mesh size. In the limit of small mesh size, the number of peaks corresponds to the number of spikes. Our comparison shows that a resolution of about 2 cm and 10 ms gives consistency between the density plot and estimated probability density distribution. It is reasonable to consider that such a scale for coarse graining of the firing density is necessarily derived from the experimental error such as tracking error (i.e., a few cm).
DISCUSSION
Complexity of the phase precession phenomenon
In hippocampal CA1 pyramidal cells, the distribution of spike phase versus position tends to be better described by a mixture of two normal density functions, N1 and N2, than by a single linear function. N1 and N2 are identified consistently among place fields of different cells in different experimental sessions. As shown above, N1 and N2 are discriminated from each other with respect to their relative position in the place field and relative phase (see, e.g., Fig. 6) as follows.
 1
 The center of N1 precedes the center of N2 along the direction of the rat's motion, but the two distributions overlap considerably in this axis.
 2
 The distributions of N1 and N2 along the theta phase axis are antiphase to one another and are much less overlapped than in the position axis. Strictly speaking, because of this antiphase relationship, it is not possible to define one component as “early” or “late;” however, according to previous usage, N2 would be considered as early in the phase cycle on the basis of the linear model in which there is a progressive retardation in phase through approximately one full cycle as the rat traverses the place field.
 3
 The range of variation of the parameters of N2 is usually larger than that of N1, which results in a variety shapes of the composite estimated spike density function in the positionphase plane. Sometimes N2 appears to exhibit a reverse phase shift or bimodality in preferred phase at a given position (as is seen in Fig. 3A); however, the data are presently insufficient to consider the possibility of such added complexity in any meaningful way.
The observation that two components provide a good fit to the data cannot alone be taken as evidence for the existence of two physiologically separate mechanisms, because a single Gaussian will never describe well the generally curved or more complex firing phase versus position profile (e.g., Fig. 1). Two Gaussians with free position, variance, and covariance parameters will almost certainly account for significantly more variance, even taking into account the additional degrees of freedom. It is important to consider, however, that the combination of two Gaussians will produce either an inflection or double maxima in the spike density distribution function. As shown in Fig. 3, when the N2 component is large in its ratio and compact enough to produce two clear peaks in the estimated spike density function, the actual spike density functions do in fact exhibit two maxima. Thus two Gaussian components reflect the shape of the actual distribution better than a unimodal but nonlinear function. The existence of two maxima in the phase versus position plots suggests that they might be observable when the average population activity is plotted as a function of phase. This effect can be clearly seen in Fig. 3B of Skaggs et al. (1996). In addition, the double maxima appear along theta phase in a given position where N1 and N2 significantly overlap with each other, as is seen clearly in the example of Fig. 3 A and in the collective distribution shown in the middle of Fig. 7.
The objection might be raised that changes in the rat's running velocity as he approached or departed from reward sites or corners might account for the general nonlinear shape of the phaseposition distributions if the phase shift were at least in part a function of time. While effects of this nature cannot be ruled out in a few isolated cases, net accelerations and decelerations are, on the average, equal and opposite, and place field distributions are essentially uniform on the track. Any effect of a systematic velocity change at specific locations would cancel out. Also, in data set 2 the rats ran on a circular platform with a randomly placed food reward, so there was no systematic trend in any place field. This data set exhibited the same bimodal effects as seen in data set 1. Furthermore, effects of velocity on phase precession have recently been described in detail by Ekstrom et al. (2001) in a study involving the effects of theNmethyldaspartate (NMDA) antagonist 1,2cyclopentenophenanthrene (CPP). CPP results in a substantial (25–30%) increase in velocity of the animal, and blocks experiencedependent changes in place field size; yet it has no effect on the general form of the phase versus position plots, which exhibit the same tendency toward bimodal distributions.
If the two components are indeed mechanistically distinct, it is likely that they are functionally different as well. We refer to the component N1 as the phase shift component, because, in this case, there is a clear phase retardation as a function of position. This component behaves as described initially by O'Keefe and Recce (1993), except that it spans only the early part of the place field and late part of the theta cycle. In contrast, component N2, which is centered in the latter half of the place field, exhibits no significant phase shift with position in average. One consequence of these properties is that firing phase can give detailed information about spatial location only within the early part of the place field.
Theoretical implications and possible mechanisms
Skaggs et al. (1996) pointed out that one consequence of the phase precession phenomenon is that, if the rat's position were reconstructed from the CA1 population activity, as done by Wilson and McNaughton (1993), but at a much finer temporal resolution, then “during each theta cycle, the reconstructed location will shift forward along the rat's trajectory for a net distance approximately the size of one average place field” (see alsoBurgess et al. 1994). The improvement of the spatial resolution by phase coding is shown to be significant (Jensen and Lisman 2000). A possible interpretation of this effect is that the phase precession reflects, in a sense, a prediction of the rat's location over the next 1–2 s. This concept has been accounted for by network models (e.g., Tsodyks et al. 1996;Wallenstein and Hasselmo 1997) in which the external input to the network was gated by the theta cycle, and hence sets the state of the network at the phase of the theta cycle at which net activity begins to rise. Asymmetry in the connections then causes the internal dynamics of the network to shift the representation further along the route during the rest of the theta cycle. This cyclic “lookahead” process was reconstructed from ensemble recordings of CA1 pyramidal cells by Samsonovich and McNaughton (1997).
The idea that the initial firing at the beginning of a place field represents a prediction of later input is supported by two further observations. First, during repeated traversals of a unidirectional path, place fields expand with experience in a direction opposite to the rat's direction of motion (Mehta et al. 1997). This effect is predicted by theoretical models of sequence learning in which potentiation of synapses from elements representing earlier parts of a sequence onto later ones leads to subsequent predictive activation of the later elements (Abbott and Blum 1996;August and Levy 1999; Jensen and Lisman 1996b; Kleinfeld and Sompolinski 1988;Levy 1989, 1996; Wallenstein and Hasselmo 1997). The range of phase precession is not affected by this expansion (Ekstrom et al. 2001; Shen et al. 1997), with the result that the slope of the phase change with position is reduced. This effect is predicted by theTsodyks et al. (1996) model as a consequence of an increase in the asymmetry parameter in the synaptic matrix. Second,Skaggs et al. (1996) showed that the spatial information conveyed by a single spike from a CA1 cell is greater during the portion of the activity identified here by N2, particularly in an open field, where place fields are not directionally dependent and a given location can be approached from more than one direction. Thus N2 associated firing appears to occur within a more restricted region of space, whereas N1 related firing might occur on any path to this region. This model can be considered to belong to a class of models that are prospective in nature, in the sense that early activity predicts subsequent input. Although it is possible that N2 is generated in CA3 or in CA1 itself by some local process, within the prospective theoretical framework, the N2 component could reflect external input (i.e., activity from other brain regions) whose availability is gated over the “early” half of the theta rhythm. The N1 component could reflect the predictive activity. If this were the case, then the observation that the N2 component is greatly attenuated or absent in DG cells implies that the external input actually arrives at subsequent processing stages. Given that CA3 and DG share a common input from layer II of the entorhinal cortex (EC), the most likely route for the external input to CA1 would be over the direct pathway from layer III of the entorhinal cortex to the CA1 pyramidal cells. This hypothesis is supported by two observations. First, Mizumori et al. (1989) showed that injection of lidocaine into the medial septum can virtually abolish both the theta rhythm and activity of CA3 pyramidal cells, while leaving placespecific firing in CA1 qualitatively unchanged. Second, EI Moser and MB Moser (personal communication) have shown that longitudinal knife cuts along the CA3CA1 border in the dorsal hippocampus severely disrupt place learning, but again leave CA1 placeselective activity qualitatively unchanged.
The implication of the foregoing theoretical interpretation is that the predictive activity (and hence the phase precession itself) should take place within or before the superficial layer of the EC and should be propagated to DG and CA3, and thence to CA1. In agreement with this idea, the phase shift component in DG is advanced in phase by 0.2 theta cycle relative to CA1. This time difference is comparable with the time for signal transmission in the feed forward pathway of the hippocampal circuit. The prospective model in general depends on asymmetry in the connections among cell assemblies corresponding to sequential experiences. Presumably these must be developed by experiencing the sequence. Hence, in the prospective model, phase precession is a consequence of sequence learning, not a mechanism to facilitate sequence learning. A major experimental challenge to this type of model would be the presence of phase precession in the activity of hippocampal pyramidal cells on the first pass through the place field in a novel environment. Preliminary observations suggest that this may be the case (Rosenzweig et al. 2000).
A second general class of theoretical models that are based on oscillator dynamics (Bose et al. 2000; Kamondi et al. 1998; O'Keefe and Recce 1993;Yamaguchi and McNaughton 1998), could be considered as retrospective in nature. These models assume the existence of oscillators with difference frequencies. O'Keefe and Recce (1993) and Bose et al. (2000) assume, in the literal sense of “phase precession,” a phase shift between slow and fast oscillators, while the former assumes a simple summation and the latter assumes transient dynamics after a trigger signal. Phase locking between nonlinear oscillators with different native frequency was proposed by Kamondi et al. (1998) based on in vitro experiments on CA1 pyramidal cells with externally applied depolarization current. Yamaguchi and McNaughton (1998)proposed a model of phase locking within a population of neurons assumed to be in entorhinal cortex or earlier. If the intrinsic frequency of a given cell increases over time following the effective stimulus, then the neuron's firing will shift to progressively earlier phases as a result of its interaction with the mean field oscillation. It is interesting that phase locking models show a phase shift ≤180°, which is consistent with the observed behavior of N1.
In the retrospective model, precession is not a consequence of learning, but may be seen as a mechanism to facilitate learning of sequences as suggested by Skaggs et al. (1996). One problem with storing sequences by sequential association is that the sequence must change rapidly on a time scale determined by the time constant of the association mechanism (e.g., t < 100 ms for NMDA receptor–dependent LTP) (Bi and Poo 1998;Levy and Steward 1983). As pointed out byMcNaughton and Morris (1986) and others, if the sequence changes too slowly, then the network performs autoassociation, and retrieval gets trapped at a fixed point. In the rat, a typical place field occupies a space of about 20–30 cm (depending somewhat on how it is measured and on behavioral and environmental factors). On linear tracks, rats typically run about 20 cm/s. In 100 ms, therefore the rat covers <1/10 of the place field, which means that the population ensemble patterns will be very highly correlated over this time window and sequence retrieval would be difficult at best. Phase precession results in overlapping segments of the recently experienced sequence being replicated repeatedly in compressed form within each theta cycle. The compression has the consequence that the patterns at the beginning and end of the cycle are essentially orthogonal. Thus a retrospective phase precession phenomenon could be seen as a device that facilitates the asymmetric strengthening of connections that may be required for sequence learning. One experimental observation in favor of retrospective encoding is the “momentum” effect described by Redish et al. (2000) in which place cells in rats running linear tracks sometimes continue to fire if the animal makes an unusual turn in the middle of the track, so that it is moving in the opposite direction. On linear tracks, firing in the opposite directions at the same location is typically uncorrelated. It should be noted that the potential generation of asymmetry is restricted to N1. In N2, the phase relations are essentially random over the range of about onehalf cycle. N2 might thus add an element of symmetry to the connections among sequentially activated neuronal populations. It is possible that N2 reflects a dominant influence of extrinsic connections during this portion of the theta cycle.
The fact that synaptic plasticity in the hippocampus may depend rather strongly on both theta phase and relative spike timing (Bi and Poo 1998; Chattarji et al. 1989; Holscher et al. 1997; Huerta and Lisman 1995;Pavlides et al. 1988; Stanton and Sejnowski 1989) would also need to be considered in any functional interpretation of the phase precession dynamics.
It should also be reemphasized that the improvement afforded by two Gaussian components in the explained variance of the spike distribution in the positionphase plane does not necessarily imply that there are two mechanistically distinct physiological processes at play. For example, the general shape of the distribution might be explainable by a simple acceleration of the sequence readout during each theta cycle. Such acceleration would result in a more rapid phase shift during the latter portion of the place field. Indeed, evidence for such acceleration can be seen in either the simulation data of Wallenstein and Hasselmo (their Fig. 7A) or the phaselocking model byYamaguchi and McNaughton (1998, their Fig. 4B).
Conclusion
In summary, description by two normal components was evaluated to quantify theta phase precession of hippocampal place cells in freely running rats. It was found that a combination of two normal density functions accounts for more of the observed variance of the actual spike density distributions than previous linear models. The two components differed in their properties. One (N1) component is localized in the initial portion of the place field and over the last half theta cycle. This component has a high negative correlation, and thus could possibly encode location in terms of theta phase. In a portion of units, a major second component (N2) occurs later in the place field, at the early half cycle and has no significant correlation on average. The phase shift component (N1) in DG is advanced by 0.2 phase units relative to that at CA1. The second component (N2) appears strongly in some but not all units in CA1 but is greatly reduced or absent in DG. These observations suggest that the phase shift component is generated at an early stage of the hippocampal circuit and that the second component is generated in CA3 or CA1, or is conveyed to CA1 via the direct entorhinal inputs. These findings were discussed in terms of their implications for both retrospective and prospective models for phase precession and its possible relation to sequence learning; however, it should be emphasized that the prospective and retrospective models are by no means mutually exclusive. It is possible that in rats running a recently experienced route, the N1 component contains activity related to both kinds of effect (see discussionabove about place field expansion). Thus further experimental work will be required to evaluate the validity of description by two components. The current analysis does show that the phase precession phenomenon may be considerably more complex than can be accounted for by any simple linear model. We believe that the further elucidation of this complexity will have important implications for understanding the neural basis of learning and memory for temporal sequences.
Acknowledgments
The authors acknowledge Drs. Carol Barnes and W. Skaggs for helpful support for the availability of the experimental data, S. Taro for cooperation on the computer programming, and Dr. Masami Tatsuno for critical reading and comments.
Footnotes

Address for reprint requests: Y. Yamaguchi, RIKEN Brain Science Institute (BSI), 21 Hirosawa, Wako, Saitama 3510198, Japan (Email:yokoy{at}brain.riken.go.jp).
 Copyright © 2002 The American Physiological Society