INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

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 7-Hz 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 large-scale 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 long-term 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.

View larger version (26K): [in this window] [in a new window]   Fig. 1. Examples of the relations between the theta phase  of spike discharge and the rat's position p for 3 individual place cells. Each dot represents 1 recorded spike. The direction of motion is from left to right. Boxes and arrows below the plots in this figure, and others mark food locations and the corners of the triangle, respectively (replotted from Skaggs et al. 1996).

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 firing-phase 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 two-component 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

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

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 and 3, cells were recorded from simultaneously in CA1 and the dentate gyrus (DG). For position analysis, light-emitting 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 by Skaggs 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); and 3) 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

(1)

where represents the theta phase (mod 1) and t is the time stamp of the ith spike of unit S and t^r_n(i)1 and t^r_n(i) are the time stamps of the n  1 and n theta peaks surrounding the spike.

The rat position p, where the ith spike of the sth unit occurred, was obtained by transformation of the two-dimensional position into a one-dimensional coordinate p on an idealized running path for each task. In the triangular maze task, the idealized running path was constructed by straight lines and circular arcs that were smoothly connected with one another so that the tangential line of the circular arc coincided with the next line. In the circular maze task, a single circle was used for the idealized path. For the linear track, the projection onto the center of the track served as the spatial coordinate. These idealizations well approximate the actual paths of the rat averaged over all laps for a given session. Because the motion of the rat was always counterclockwise along the rectangular track and circle, increasing values of p corresponds to forward motion of the rat. For the linear track, the origin of the forward and return journeys were considered as the zero reference, and the position coordinate was defined in terms of positive distance from the origin of the current journey. It is well-documented that on-linear tracks, place fields in the forward, and return directions are uncorrelated (Markus et al. 1994; McNaughton et al. 1983; Muller et al. 1994). The plot of unit firing in the p- plane is referred to below as a p- plot. For both illustrative and analysis purposes, the  axis in some cases is plotted over several cycles.

Estimation of the probability density function by means of the expectation maximization (EM) algorithm

For the purpose of analysis, the spike distribution in the p- 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 best-fitting 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

(2a)

with

(2b)

and

(2c)

In the present application, x and y correspond with the position p and theta phase  at which the spike event occurs. The parameters, µ_x and µ_y represent the mean values with respect to x and y coordinate, and , , and _xy denote the variance of x, y, and their covariance, respectively. The quantity r is a correlation coefficient and ranges between 1 and 1.

When the probability density function is assumed to consist of N normal density functions, the function P(x, y) becomes

(3a)

with

(3b)

where suffix i denotes the ith normal density function (i = 1 ~ N). The factor _i denotes the relative weight of the ith normal density function.

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 two-dimensional 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, and L normal density functions are assumed in a cycle, the EM algorithm is applied with n = LM in Eq. 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.).

View larger version (41K): [in this window] [in a new window]   Fig. 2. An example of unit activity distribution in the p- plane and an introduction to the p- coordinate system used in the expectation maximization (EM) algorithm. The actual spike distribution from 1 unit is shown in A and the corresponding density plot is shown in B. The EM algorithm with L = 2 was applied to these data plotted over 5 theta cycles in the p- plane, and the results are shown in C and D. C: the contour lines of individual normal density functions at the level of 1/e. D: the contour plot of the probability density function P given by 10 normal density functions. The middle cycle of the  coordinate [2 + 0, 3 + 0] is used for further statistical evaluation below (0 = 0.25). The horizontal bar in the bottom right side represents 10 cm.

The new variable  representing M periods of theta phase is introduced below.  is not a modulus but a real number defined in the interval [₀, ₀ + M]. The value ₀ 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., [₀ + 1, ₀ + 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 and 2 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.25. The result of applying the EM algorithm, with L = 2, to this unit is shown in Fig. 2, C and D. Figure 2C shows the 10 individual normal density functions, each of which is denoted by a contour line at the level of 1/e. Figure 2D 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} for L functions. Thus the estimated parameters with a given value of L are obtained for each unit as follows

(4a)

with

(4b)

Now we have the results of the EM algorithm as {_L}, which determines the probability density functions for individual units. The correlation coefficient termed r_{p j} is given by Eq. 2c with _L. Thus the obtained probability density function is written as P(p,   _L).

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

(5a)

where (p_k, _k) denotes the value of (p, ) for an observed spike with a tentative number k in the set of all spikes in each unit. P(p,  *) denotes the conditional probability. Either P(p_k, _k L) or P(p_k, _k *) represents the value of the probability density for the kth spike. In the present case, the conditional probability is given by the function with the maximum number of L, L = 3. That is

(5b)

To evaluate the collective properties of the probability density functions in the set of units, two functions are used. One is the probability density function with parameters obtained by averaging the parameters of µ_p1 or µ_p2 over the set of units. The other is the average of probability density functions over the set of units, after aligning the distributions on either µ_p1 or µ_p2 of the individual units.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

General results

Several examples of results obtained in data set 1 by the EM algorithm are shown in Fig. 3. The results with L = 1, 2, and 3 are shown for the sake of comparison. In Fig. 3A, 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 "S-shaped" 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.

View larger version (23K): [in this window] [in a new window]   Fig. 3. Three examples of the results obtained by the EM algorithm are shown. A, top: actual spike distribution (left) and corresponding density plot (right). The next 3 rows of plots illustrate the results obtained with L = 1, 2, and 3, respectively. The contours at 1/e of the peak of individual normal density functions are shown on the left, and the contour plots of the density distribution functions are shown on the right. Note that the density functions for both L = 2 and L = 3 contain 2 peaks. Horizontal scale bar shows 10 cm. B: data from a different unit plotted as in A. Note that the density function again shows 2 peaks for both L = 2 and L = 3. C: an example of a unit with a probability density function consisting of a single peak for L = 1, 2, and 3.

Figure 3B 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 and L =3.

Figure 3C 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 or B.

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 for L = 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 for L = 2 in a half of units. The result for L = 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 ₁ and ₂ is larger than that between ₂ and ₃ 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.

View larger version (10K): [in this window] [in a new window]   Fig. 4. The likelihood parameters L with L = 1, 2, and 3 are plotted for individual units in data sets 1 and 2. A: CA1 units (n = 29) in data set 1. B: dentate gyrus (DG) units (n = 4) in data set 1. C: CA1 units (n = 15) in data set 2. Note that there is substantial variation in the improvement in fit obtained in going from L = 1 to L = 2 but that there is little further improvement at L = 3.

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 and its SD are 0.49 ± 0.12 in data set 1 and 0.51 ± 0.16 in data set 2. That is, the relative weights are almost equal to each other.

Figure 6 summarizes the results of CA1 units in data sets 1 and 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 in data set 1 and 12.98 ± 4.30 cm in data set 2. In data set 1 shown in Fig. 5A, the distributions of µ₁ and µ₂ 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 anti-phase to each other. The anti-phase relation between N1 and N2 is similarly observed in data set 2 as shown in Fig. 5B.

View larger version (18K): [in this window] [in a new window]   Fig. 5. Results of analysis of CA1 units in data set 1 (n = 29) and data set 2 (n = 15) obtained with L = 2. A: relative location of 2 centers (N1, N2) in data set 1. The centers of N1 and N2 are plotted as open and closed circles, respectively. For the sake of comparison, the origin of the position coordinates for each unit are aligned on the center of N1. B: relative location of 2 centers in a session of data set 2. C: correlation coefficients of N1 and N2 in data set 1. D: correlation coefficients of N1 and N2 in data set 2. E: gradient of principal axis of N1 of all CA1 units in data set 1. F: gradient of principal axis of N1 of all CA1 units in the session of data set 2.

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 anti-phase 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.

View larger version (35K): [in this window] [in a new window]   Fig. 6. The distributions of N1 and N2 along the p coordinate (A) and along the  coordinate (B) are plotted for every CA1 unit in data set 1 (n = 29). In A, the closed circles show the center of N1 (aligned as µpi  µp1), and the solid line shows the width of N1, p1. The open squares show the relative center of N2, and the shaded line represents the width of N2, p2. In B, the center and width of N1 and N2 are shown by the same symbols as in A. The distributions of N1 and N2 for every CA1 unit in data set 2 (n = 15) are plotted in C and D, respectively.

The relation between position and theta phase in N1 and N2 are summarized in the lower part of Fig. 5. The correlation coefficients r_pj 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 {₂} 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.

View larger version (39K): [in this window] [in a new window]   Fig. 7. Probability density distributions obtained by averaging the estimated parameters of individual units in data set 1. A: average of all CA1 units. B: average of CA1 units with double peaks. C: average of CA1 units with a single peak.

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, anti-phase 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. 7A 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.

View larger version (33K): [in this window] [in a new window]   Fig. 8. The probability density function obtained by the average of estimated parameters in data set 2. A: average of all CA1 units (n = 15). B: average of CA1 units with double peaks (n = 13). C: average of CA1 units with a single peak (n = 2).

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, A and 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 coefficient r 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.

View larger version (18K): [in this window] [in a new window]   Fig. 9. Results of analysis of DG units in data sets 1 and 3. A: relative location of 2 centers (N1, N2) in data set 1 (n = 4). B: relative location of 2 centers in a session of data set 3 (n = 9: since 0 are different with each experimental session in data set 3, only 1 session that shows most DG units is analyzed). C: correlation coefficients of N1 and N2 in data set 1 (n = 4). D: correlation coefficients of N1 and N2 in data set 3 (n = 20). E: gradient of principal axis of N1 of all DG units in data set 1 (n = 4). F: gradient of principal axis of N1 of all DG units in the session of data set 3 [n = 19, 1 unit had a low correlation coefficient (<0.1) and is omitted from the figure].

View larger version (16K): [in this window] [in a new window]   Fig. 10. Results of analysis of CA1 units (n = 10) in the same session as that of data set 3 in Fig. 9. A: the correlation coefficients of N1 and N2. B: histogram of the gradient of the principal axis of N1.

Figure 11 shows the probability density functions obtained as an average of the parameters {₂}, for DG units. The average in data set 1 (Fig. 11A) 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.

View larger version (79K): [in this window] [in a new window]   Fig. 11. The probability density functions of DG units obtained from the average of estimated parameters. A: all DG units (n = 4) in data set 1. B: all DG units (n = 9) in data set 3. C: DG units with double peaks in the session of data set 3 (n = 4). D: DG units with a single peak in data set 3 (n = 5).

Since one-half 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. 11B, over units with double peaks in Fig. 11C and over units with a single peak in Fig. 11D. It is seen that all these results show very similar distributions with a single peak. The second peak of individual units disappeared in Fig. 11C 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 11A. 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. 12B) 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.

View larger version (87K): [in this window] [in a new window]   Fig. 12. Summation of the probability density functions over a set of units in data set 1 is shown by contour lines. In the summation, the data are aligned on the center of N1. A: all CA1 units (n = 29). B: CA1 units with double peaks (n = 23). C: CA1 units with a single peak (n = 6). D: all DG units (n = 4).

View larger version (91K): [in this window] [in a new window]   Fig. 13. Summation of the probability density functions over a set of units in data set 1 is shown by contour lines. In the summation, the data are aligned using the center of N2. A-D are as in Fig. 12.

For the sake of comparison, Fig. 14 shows 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).

View larger version (99K): [in this window] [in a new window]   Fig. 14. Density plots of spikes in the original experimental data of CA1 and DG units in data set 1. Spike density is averaged using a mesh of 1.56 [cm] × 0.056 [theta cycle] and represented by 15 linearly spaced contour lines. A: the spike densities are summarized in a group of CA1 units with a single peak (n = 6), using N1 center as the origin of position coordinate. B: the spike densities are summarized in a group of CA1 units with double peaks (n = 23), using N1 center as the origin of position coordinate. C: the same data as in A are shown using N2 center as the origin of position coordinate. D: the same data as in B are shown using N2 center as the origin of position coordinate. E: data from all DG units plotted using N1 center as the origin of position coordinate (n = 4). F: same data as in E plotted using N2 center as the origin of position coordinate.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

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 anti-phase to one another and are much less overlapped than in the position axis. Strictly speaking, because of this anti-phase 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 position-phase 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. 3A 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 phase-position 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 the N-methyl-D-aspartate (NMDA) antagonist 1,2-cyclopentenophenanthrene (CPP). CPP results in a substantial (25-30%) increase in velocity of the animal, and blocks experience-dependent 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 also Burgess 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 "look-ahead" 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 the Tsodyks 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 place-specific firing in CA1 qualitatively unchanged. Second, EI Moser and MB Moser (personal communication) have shown that longitudinal knife cuts along the CA3-CA1 border in the dorsal hippocampus severely disrupt place learning, but again leave CA1 place-selective 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