## Abstract

Spatial information about the environment is encoded by the activity of place and grid cells in the hippocampal formation. As an animal traverses a cell's firing field, action potentials progressively shift to earlier phases of the theta oscillation (6–10 Hz). This “phase precession” is observed also in the prefrontal cortex and the ventral striatum, but mechanisms for its generation are unknown. However, once phase precession exists in one region, it might also propagate to downstream regions. Using a computational model, we analyze such inheritance of phase precession, for example, from the entorhinal cortex to CA1 and from CA3 to CA1. We find that distinctive subthreshold and suprathreshold features of the membrane potential of CA1 pyramidal cells (Harvey et al., 2009; Mizuseki et al., 2012; Royer et al., 2012) can be explained by inheritance and that excitatory input is essential. The model explains how inhibition modulates the slope and range of phase precession and provides two main testable predictions. First, theta-modulated inhibitory input to a CA1 pyramidal cell is not necessary for phase precession. Second, theta-modulated inhibitory input on its own generates membrane potential peaks that are in phase with peaks of the extracellular field. Furthermore, we suggest that the spatial distribution of field centers of a population of phase-precessing input cells determines, not only the place selectivity, but also the characteristics of phase precession of the targeted output cell. The inheritance model thus can explain why phase precession is observed throughout the hippocampal formation and other areas of the brain.

## Introduction

Place-specific firing and local field potential (LFP) oscillations in the rodent hippocampal formation could be important to understanding navigation and episodic memory (Buzsáki, 2005). Within a firing field, action potentials exhibit phase precession; that is, the firing phases relative to the LFP theta decrease as a function of position (O'Keefe and Recce, 1993; Hafting et al., 2008). Phase precession can encode spatial information (Jensen and Lisman, 2000; Reifenstein et al., 2012a) and can temporally compress behavioral sequences to match the time scales of neuronal plasticity (Skaggs et al., 1996; Bi and Poo, 1998). Mechanisms generating phase precession, however, remain unresolved (Maurer and McNaughton, 2007; Burgess and O'Keefe, 2011; Eggink et al., 2014). Network models of phase precession emphasize the role of recurrent inputs from cells with neighboring place fields (Jensen and Lisman, 1996; Tsodyks et al., 1996; Wallenstein and Hasselmo, 1997; Baker and Olds, 2007; Navratilova et al., 2012), coupled excitatory and inhibitory neurons (Bose et al., 2000; Castro and Aguiar, 2012; Cutsuridis and Hasselmo, 2012), and a summed-population effect (Geisler et al., 2010; Thurley et al., 2013). Alternatively, cellular models suggest that phase precession results from an interaction between a somatic and a dendritic signal within a cell (Kamondi et al., 1998; Magee, 2001; Harris et al., 2002; Mehta et al., 2002; Leung, 2011), a superposition of two oscillations with different frequencies (O'Keefe and Recce, 1993; Lengyel et al., 2003; O'Keefe and Burgess, 2005), synaptic facilitation (Thurley et al., 2008), persistent firing (Hasselmo, 2008), or excitatory inputs from two regions (Chance, 2012). The possibility of feedforward inheritance, however, has not been studied thoroughly.

Phase precession is observed, not only throughout the hippocampus (for example, see Skaggs et al., 1996), but also in the medial entorhinal cortex (MEC; Hafting et al., 2008), the prefrontal cortex (Jones and Wilson, 2005), the subiculum (Kim et al., 2012), and the ventral striatum (van der Meer and Redish, 2011). Interestingly, principal cells in these regions receive projections from regions that also show phase precession. Can phase precession observed in a cell be simply inherited from a population of cells upstream?

In this computational study, we investigate how phase precession can be propagated from one region to another. We first analyze how phase precession in CA3 can lead to characteristic subthreshold (Harvey et al., 2009) and suprathreshold (for example, see Mizuseki et al., 2012) features of phase precession in a CA1 place cell. We also outline how inhibition modulates phase precession in CA1 (Royer et al., 2012). Next, in a more general feedforward network, we focus on spatial distributions of field centers of excitatory phase-precessing input cells, including grid cells in the MEC. We determine the resulting place selectivity and the phase precession characteristics of a target cell and explain, for example, why phase precession in hippocampal interneurons is weaker than in principal neurons (Maurer et al., 2006a; Ego-Stengel and Wilson, 2007; Geisler et al., 2007). The proposed inheritance mechanism may thus account for phase precession in various cell types in the hippocampus proper, subiculum, prefrontal cortex, and ventral striatum.

## Materials and Methods

##### Model of the CA1 membrane potential.

We model the membrane potential *V*_{CA1}(*t*) at the soma of a CA1 pyramidal cell by taking into account an inhibitory and an excitatory component. Inhibitory input is characterized by an ongoing oscillation of the membrane voltage of a CA1 pyramidal cell. The oscillation is coherent with the LFP (Kamondi et al., 1998; Bland et al., 2005) and is independent of place field activity (Harvey et al., 2009; Epsztein et al., 2011). We model this intracellular theta oscillation with a sinusoidal function as follows:
where *B* is the amplitude of the oscillation, ω_{θ} is its angular frequency (ω_{θ} = 2π*f*_{θ} and *f*_{θ} is the frequency of the theta LFP), and ϕ_{θ} is the phase of the intracellular oscillation peak with respect to the theta LFP peak. Here, we use the convention that 0° corresponds to the peak of the theta LFP recorded in the CA1 pyramidal cell layer. We note that the DC level of the inhibitory contribution *V*_{θ} (*t*) is below zero; that is, the oscillation is always hyperpolarizing.

Phase-precessing excitatory input from CA3 to CA1 is characterized by the dynamics of the CA3 population firing rate and the shape of EPSPs; input from the entorhinal cortex is analyzed at a later stage. To model the excitatory component from CA3, we assume that the animal is traversing a place field with constant speed so that time and space are equivalent up to a constant conversion factor. The place fields of the cells of the CA3 population are assumed to be identical and completely overlapping. This assumption will be relaxed later.

The time course of the firing rate of a CA3 cell is described by two parts. The first part reflects the place field: the average firing rate increases when the animal enters the place field, achieves its maximum at the center of the place field, and decreases when the animal exits the place field. We describe this ramp-like behavior of the firing rate with a Gaussian function *x*(*t*) around the center *t _{c}* of the place field as follows:
where σ determines the place-field width (∼3σ). The second part of the firing rate model is an oscillatory modulation that accounts for CA3 phase precession. The oscillatory modulation is determined as follows:
and is characterized by the angular frequency ω

_{λ}, where ω

_{λ}= 2π

*f*

_{λ}and

*f*

_{λ}is the frequency of the oscillation, the phase offset ϕ

_{λ}with respect to the reference theta LFP, and the modulation depth

*C*. Multiplying this oscillatory modulation with the ramp

*x*(

*t*) in Equation 2, we find the firing-rate model of the input of one CA3 cell to a CA1 cell as follows: where λ

_{0}is the average firing rate in the center of the field. The combined firing rate from a population of

*N*identical CA3 cells is then

*N*λ(

*t*).

To ensure that the oscillating, place-field like excitatory input λ(*t*) to a CA1 pyramidal cell shows phase precession, we assume that *f*_{λ} is greater than the theta frequency *f*_{θ} = 8 Hz (Geisler et al., 2010); furthermore, we require several oscillation cycles within a place field, that is, 3σ ≫ 1/*f*_{θ}; typically σ ∼ 0.3 s corresponding to a place-field size of ∼1 s (Fig. 1*A*).

The rate λ(*t*) serves as a basis to generate the spiking input to a CA1 pyramidal cell. To model the variability of spiking input, we describe the activity of a population of *N* CA3 cells by an inhomogeneous Poisson process with time-dependent firing rate *N*λ(*t*). We represent a spike at time *t _{f}* as a delta function δ(

*t*−

*t*) centered at the time of firing. Combining the spikes from all

_{f}*N*input cells, we obtain the total neural response function as follows: where

*t*is the

_{f}*f*th spike time. An example is shown in Figure 1

*B*. For the analysis and the simulations, it is important at what times neurons in the input population were active, but it is not important which particular CA3 neuron actually fired. Therefore, we do not keep track of the CA3 neuron that fired.

To describe how an input spike affects the membrane potential of a CA1 pyramidal cell, we assume that each spike elicits an EPSP that is identical for all spikes and all input neurons. For simplicity, we model EPSPs by an α function as follows:
where ϵ_{max} is the maximum amplitude obtained at *t* = τ and τ determines both the rise time and decay time of the EPSP. We restrict the domain of ϵ(*t*) to positive time values, that is, *t* > 0, and we assume that ϵ(*t*) = 0 for t ≤ 0. The CA3 contribution to the CA1 excitation *v*_{input}(*t*) for a Poisson spiking model is then calculated as follows:
where * represents the operation of convolution. An example simulation is shown in Figure 1*C*.

To model the time course of the subthreshold membrane potential of a CA1 pyramidal cell during a place-field traversal, we add the inhibitory component in Equation 1 and the excitatory component in Equation 8. The total subthreshold membrane potential is calculated as follows:
where *V*_{rest} is the resting membrane potential, which is −70 mV in the simulations. The example simulation in Figure 1*E* shows that outside the place field, that is, where excitatory input is negligible, the membrane potential oscillates at theta frequency *f*_{θ} with amplitude *B* and peaks at phase ϕ_{θ} = 0° with respect to the LFP (i.e., peaks concurrently with the LFP peaks).

Within the place field, inhibitory and excitatory input interfere. The phases (again with respect to the LFP) of local maxima of the membrane potential are plotted against time in Figure 1*F*.

##### Analytical expressions for the mean-field subthreshold membrane potential.

To mathematically describe properties of the model of the membrane potential of a CA1 pyramidal cell, we first derive expressions for the mean field; that is, the trial-averaged membrane potential. The noise due to the variable and discrete spiking activity of the CA3 neurons is characterized at a later stage.

The mean-field expression for the excitatory input in Equation 8 is as follows:
where 〈.〉 indicates trial averaging. Note that we have used an uppercase *V* to denote a mean-field voltage. To proceed, we consider the kernel ϵ(*t*) as a filter ϵ̄(ω) := *F*(ϵ(*t*))(ω) in the frequency domain where *F* denotes the Fourier transform. The filter ϵ̄ has magnitude |ϵ̄(ω)| = (eϵ_{max}τ)/(1 + ω^{2}τ^{2}) and phase ϕ(ω) = arg(ϵ̄(ω)) = 2tan^{−1}(ωτ) where arg denotes the operation that yields the angle of a complex number. Using Equations 4 and 6, we can approximate *V*_{input}(*t*) for ω = ω_{λ} ≫ 1/σ as follows:
where τ_{ϕ} = −(∂ϕ(ω))/(∂ω) is the phase delay and τ* _{g}* = −ϕ(ω)/ω is the group delay (Haykin and Moher, 2009). For the time scales involved in this integration problem (ω

_{λ}τ ≲ 1), we can approximate τ

_{ϕ}≈ τ

*≈ 1.8τ and rewrite the above expression as follows: The excitatory contribution*

_{g}*V*

_{input}(

*t*) to the CA1 membrane voltage is therefore a scaled (by e

*N*ϵ

_{max}τ), temporally shifted (by 1.8τ), and filtered (by 1/(1 + ω

^{2}

_{λ}τ

^{2})) version of the CA3 firing rate λ(

*t*) in Equation 4. The total mean-field membrane potential, including the ongoing theta oscillation and the resting membrane potential, is as follows: To further characterize the mean-field membrane potential

*V*

_{CA1}and to be able to compare its properties to experimentally accessible quantities, we assume that, in the center of the place field (where

*x*(

*t*) ≈ 1), the oscillation of the membrane potential is mainly due to the excitatory input from CA3, which means that the oscillation amplitude due to the ongoing inhibitory input can be neglected. In other words, we assume that e

*CN*λ

_{0}ϵ

_{max}τ ≫

*B*. Using this approximation, in the center of the place field the oscillation amplitude is calculated as follows: Furthermore, the CA1 membrane potential shows a mean depolarization ramp that is calculated as follows: We can readily calculate the modulation depth

*C*of the CA1 membrane potential by taking the ratio of the oscillation amplitude Δ

_{VCA1}*V*

_{osc}and the average depolarization Δ

*V*

_{ramp}as follows: where

*C*is, as before, the modulation depth of the CA3 population activity. In other words, the modulation depth

*C*of the CA1 mean-field membrane potential oscillation is a filtered (by 1/(1 + ω

_{VCA1}^{2}

_{λ}τ

^{2})) version of the modulation depth

*C*of the CA3 firing rate.

##### Quantitative analysis of membrane-potential noise.

The membrane potential of the CA1 model cell shows fluctuations because the input is described by discrete spikes and spike times are generated by a Poisson process. Below, we perform a quantitative analysis of this variability that is termed shot noise. The properties of this noise will then be compared with the properties of the mean-field membrane potential.

One way of quantifying the effects of shot noise on the excitatory contribution *v*_{input}(*t*) is by calculating its SD as follows:
where 〈.〉 denotes, again, an average of realizations, or trials, of the Poisson process. We note that this SD of the membrane potential does not take into account the ongoing oscillations *V*_{θ}.

The second moment of the membrane potential in Equation 8 is defined as follows:
The next step in the analysis is to separate the terms dependent on the realization from those that are independent. For this, we partition the spike train into time intervals *I _{l}* = (

*t*,

_{l}*t*

_{l}_{+1}) of infinitesimal length Δ

*t*. We introduce a function Γ(

_{l}*l*) that counts the number of events (i.e., spikes) in the interval

*I*. We can thus write the second moment as follows: where we have separately considered the cases

_{l}*l*≠

*l*′ and

*l*=

*l*′. Using the definition of an inhomogeneous Poisson process (for example, see Kempter et al., 1998), we perform an average of the functions Γ over realizations to obtain 〈Γ(

*l*)Γ(

*l*′)〉 =

*N*

^{2}λ(

*t*)Δ

_{l}*t*λ(

_{l}*t*)Δ

_{l′}*t*and 〈Γ

_{l′}^{2}(

*l*)〉 =

*N*λ

*Δ*

_{l}*t*. With these averages, we can express Equation 22 as a Riemann sum that equates to two convolutions in the limit when Δ

_{l}*t*and Δ

_{l}*t*tend to zero: Finally, using Equations 12 and 23, we can calculate the variance of the membrane potential as follows: To obtain a closed expression of the variance, we approximate the function

_{l′}*t*− τ) depends on time

*t*, as does the variance 〈Δ

*v*

^{2}

_{input}〉. To obtain a time-independent estimate of the variance that characterizes the noise within the place field, we replace the time-dependent rate λ(

*t*− τ) by the average firing rate λ

_{0}in the center of the place field. A typical value of the SD δ

*v*

_{input}of the membrane potential, then, is as follows: We can now construct a “quality parameter” ρ of the membrane potential that characterizes the relative strength of the oscillation amplitude Δ

*V*

_{osc}with respect to the noise amplitude δ

*v*

_{input}as follows: The prefactor e/2, which is close to 1, was introduced to arrive at a simple expression devoid of an irrelevant numerical prefactor. The contributions to the quality parameter ρ are therefore: the number

*N*of active CA3 cells projecting to a single CA1 cell, the time constant τ of the EPSP, the modulation depth

*C*of the population rate oscillation at angular frequency ω

_{λ}, and the mean firing rate λ

_{0}in the center of the place field. Using Equations 16 and 17, we can rewrite the above expression of the quality parameter ρ in terms of experimentally measurable voltages Δ

*V*

_{ramp}, and Δ

*V*

_{osc}as follows: Furthermore, using Equation 16 for Δ

*V*

_{osc}, Equation 17 for Δ

*V*

_{ramp}, and Equation 29 for ρ, we can determine

*C*,

*N*, and ϵ

_{max}. To obtain an explicit expression for

*C*, we use the ratio of Δ

*V*

_{osc}and Δ

*V*

_{ramp}and rearrange the terms as follows: Using this result for

*C*in Equation 29, we find: Inserting this result for

*N*in Equation 17, we obtain: We are thus able to derive and predict values of the model parameters

*C*,

*N*, and ϵ

_{max}. To provide an estimate of the error of the prediction, for example of

*N*, based on known and measurable parameters as in Equation 32, we can use a rule of propagation of errors as follows: where δ

*N*is the (predicted) error of our estimate of

*N*and δΔ

*V*

_{ramp}, δΔ

*V*

_{osc}, δλ

_{0}, and δτ are SEMs obtained from previously published studies (Magee and Cook, 2000; Harvey et al., 2009; Mizuseki et al., 2009, 2012).

##### Suprathreshold model of phase precession.

To account for phase precession of spikes of CA1 cells, we used the two-compartment model that was developed by Pinsky and Rinzel (1994) for CA3 cells and adapted to the CA1 region by Kamondi et al. (1998). This particular model reproduces bursting in CA1 pyramidal neurons. The voltages *V _{s}* and

*V*are described by the following differential equations: where

_{d}*s*and

*d*are the somatic and dendritic compartments, respectively,

*C*is the membrane capacitance per unit area,

_{m}*I*is a (ohmic) leak current density,

_{L}*I*

_{Na}is a sodium current density,

*I*

_{K}is a potassium current density,

*I*

_{Nap}is a persistent sodium current density,

*I*

_{Ks}is a slow potassium current density,

*g*is the coupling conductance between the two compartments (

_{c}*s*and

*d*),

*p*is the ratio of the somatic area to the total area,

*I*is the synaptic current density applied to the soma, and

_{s}*I*is the synaptic current density applied to the dendrite. The current densities

_{d}*I*

_{Na},

*I*

_{Nap},

*I*

_{K}, and

*I*

_{Ks}, are described by the Hodgkin-Huxley formalism (Pinsky and Rinzel, 1994). The values for the intrinsic conductances are taken from Kamondi et al. (1998).

The synaptic current density *I _{d}* in our model represents the phase-precessing activity of the presynaptic CA3 cells and is modeled as a sum of EPSCs in a similar fashion as for the subthreshold model, where the CA3 input was modeled as a sum of EPSPs. As for the subthreshold model, we assume that spikes are generated via an inhomogeneous Poisson process and we convolve the respective neural response function with an EPSC kernel. A single EPSC (kernel) is represented by an α function (as in Eq. 6) as follows:
where the subscript

*I*indicates a current kernel. The convolution of the EPSC kernel ϵ

*with the neural response function*

_{I}*R*(

*t*), that is, the CA3 spike train generated via a Poisson process, is precisely the current density

*I*that is fed onto the dendrite as follows: where

_{d}*A*is the area of the dendritic compartment. The synaptic current density

_{d}*I*, which is fed onto the soma, models the effect of the ongoing theta oscillation that is coherent with the extracellular LFP as follows: where ϕ

_{s}*is a phase offset with respect to the theta LFP. Values of the model parameters used for the simulations are specified in the caption of Figure 3.*

_{I}##### Model for place-selective responses.

In previous analyses, we assumed that the place fields of the cells of the CA3 population that project to a single CA1 cell are identical. We now relax this assumption and consider the case where the place fields of the input population are spatially distributed and partially overlapping. For this, we consider a virtual rat running on a linear track of length *L* at a constant speed *v*, so that the time to complete a run is *T*_{tot} = *L*/*v*. We assume that, along a particular direction (e.g., from left to right), there are *N* active place cells with place fields that overlap on the track. Furthermore, we assume that these cells exhibit phase precession within their place fields. These *N* input cells project to a single output cell that linearly integrates the inputs to produce a somatic voltage. We characterize the behavior of this output cell in terms of place-field selectivity and phase precession.

The rate λ(*t*,*T _{i}*) of a place cell that is active on the linear track is as follows:
where

*T*is the center of the place field,

_{i}*f*

_{λ}is the frequency of the oscillation modulating each cell,

*C*is the firing-rate modulation depth, τ

*is the theta-time scale shift of the oscillation, and σ is proportional to the place field's size. Because we assume phase precession at the input, the frequency of the oscillation is greater than the frequency of the theta LFP (*

_{i}*f*

_{λ}>

*f*

_{θ}). As shown in previous studies, the variables

*T*and τ

_{i}*are correlated so that one can write τ*

_{i}*=*

_{i}*kT*, where

_{i}*k*is a constant, the compression factor (Skaggs et al., 1996; Dragoi and Buzsáki, 2006; Geisler et al., 2010).

The population firing rate λ(*t*) is obtained through a sum over the *N* cells as follows:
For N ≫ 1 and if *T _{i}*

_{+1}−

*T*≪ σ ∀

_{i}*i*, we can introduce the density

*p*of place field centers such that the population input rate is as follows: For place-field centers within the interval

*t*) analytically once we have defined the density

*p*of place-field centers at the input. Here, we consider four distributions, namely

*p*(delta),

_{D}*p*(Gaussian),

_{G}*p*(uniform), and

_{U}*p*(ramp) defined as follows: where σ

_{R}*≪*

_{d}*T*

_{tot}is the width of the Gaussian distribution. Finally, the cell's output activity, that is, its somatic voltage

*V*

_{out}(

*t*), is calculated as follows: where ϵ(

*t*) is, as before, a kernel representing an EPSP. We have simulated the results for

*V*

_{out}(

*t*) for each of the above distributions in Figure 5

*A–D*. Here, we analytically derive λ(

*t*) for the case of

*p*=

*p*. In this case, the integral for the population rate, Equation 42, becomes: where the limits at infinity are taken provided the size σ

_{G}*of the place-field center distribution*

_{d}*p*is much smaller than

_{G}*T*

_{tot}. Solving the above integral, we obtain: The above expression represents an output place field with width (proportional to) σ

*and modulated with a frequency*

_{R}*f*given as follows: Furthermore, the modulation depth

_{R}*C*

_{out}of the output cell is calculated as follows: Now, we calculate the phase-precession range for this output place field. The range Δϕ is given as follows: where we have used the fact that

*f*

_{θ}=

*f*

_{λ}(1 −

*k*) (Geisler et al., 2010). In the limit of a large place field distribution (σ

*≫ σ), we recover the results predicted for a uniform distribution*

_{d}*p*=

*p*at the input; that is, a theta oscillation at the output with no place selectivity.

_{U}##### Grid-to-place phase-precession inheritance model.

To study inheritance from grid cells to place cells in the hippocampus, we consider the case of phase precession in grid cells of the MEC, which project to pyramidal cells in the hippocampus via the perforant path. Solstad et al. (2006) modeled the grid-to-place transformation as a sum of periodic modes, corresponding to the spacings of the grids, to produce a nonperiodic place field. We extend these results and include an oscillatory modulation of the input grid fields to account for phase precession.

To model the periodic modes, we focus on the one-dimensional case such that the grid-field activity *G* is periodic along the *x*-axis as follows:
where *G*_{max} is the maximum firing rate within the field and *s* is the grid-field spacing; that is, distance between two consecutive grid maxima. By linearly combining grid fields with different spacings *s*, we obtain a target function *P*(*x*) that represents the activity of a place field and is parametrized as follows:
where *P*_{max} is the maximum firing rate and σ determines the size of the place field. We can decompose the even function *P*(*x*) into Fourier modes as follows:
where *k* = 2π/*s* is the spatial frequency and
is the Fourier transform of *P*(*x*). In general,

To relate the place-field activity *P*(*x*) to the grid-field activity *G*(*s*, *x*) in Equation 57, we express the cosine term in Equation 59 in terms of the grid functions as follows:
The expression for *P*(*x*) in Equation 59 then reads:
where *s* = 2π/*k* and d*s* = (−2π/*k*^{2})d*k*. We can express the above equation as a Riemann sum by considering *N* ≫ 1 grid fields and by discretizing *s* as follows:
where *n* ∈ {1, 2, …, *N*} and *s*_{min} and *s*_{max} are the minimum and maximum spacings, respectively. Equation 62 can be approximated (including a factor 2 to account for both positive and negative spatial frequencies) as follows:
where *A*(*s _{n}*, σ) is the synaptic weight connecting a grid cell with grid field activity

*G*(

*s*,

_{n}*x*) to the target place cell. We note that the weight

*A*(

*s*, σ) depends on both the spacing

_{n}*s*of the input grid functions and the size σ of the target place field; see Figure 6

_{n}*E*for an example of

*A*(

*s*, σ) for σ = 0.22 m,

_{n}*N*= 50,

*s*

_{min}= 0.1 m, and

*s*

_{max}= 4 m.

To account for phase precession, we introduce an oscillatory modulation *M _{n}* of the grid-cell firing rates

*G*(

*s*,

_{n}*x*) in a similar manner as we did for CA3 cells in Equation 3. The time-dependent modulation

*M*for a grid cell with spacing

_{n}*s*is as follows: where

_{n}*C*∈ (0,1] is the modulation depth, ω

*is the angular frequency, and ϕ is a phase offset. The modulated grid activities*

_{n}*G*(

_{M}*s*,

_{n}*x*) are then: where, as in previous sections of the manuscript, we assume that our virtual animal is running at a constant speed

*v*such that

*t*=

*x*/

*v*.

Our requirements for choosing ω* _{n}* and ϕ are: (1) that there is phase precession per se (Hafting et al., 2008), (2) that the phase-precession slope is inversely proportional to the size of the firing field (Brun et al., 2008), and (3) that there is phase locking of the first spike(s) within the field (Mizuseki et al., 2009). The first requirement is fulfilled if we take ω

*≳ ω*

_{n}_{θ}∀

*n*, where ω

_{θ}is the angular frequency of the theta LFP. The second requirement is equivalent to demanding that the phase-precession range Ω is constant and the same for all grid fields; that is, the range is independent of the grid spacing. For this, we must first define the size of a single grid firing field, which in our model is closely related to the spacing

*s*. A firing rate >20% of the peak rate delineates the extent of a firing field (Brun et al., 2008; Kjelstrup et al., 2008). Using this criterion and Equation 57, we can deduce that the size of a grid firing field is ≈ 0.7

_{n}*s*, and therefore 0.7

_{n}*s*/

_{n}*v*is the time spent traversing a grid field. Therefore, the phase-precession range Ω can be expressed as follows: which implies that The simplest way to model phase locking of the first spike within a grid field and to fulfill our third requirement is to ensure that the theta phase at entry ϕ

_{entry}is the same for two neighboring firing fields, which are separated by the spacing

*s*. The difference between the entry phases (with respect to theta LFP) of the two consecutive grid fields should be 2π. This condition is equivalent to the following: Let us express the phase offset ϕ in Equation 65 in terms of the experimentally measured entrance phase ϕ

_{n}_{entry}. Because of the phase-locking condition, the entry phase ϕ

_{entry}is the same for neighboring grid fields. Therefore, the relationship between the entrance phase ϕ

_{entry}and the phase offset ϕ can be determined for any firing field. Focusing on, for example, the central grid field, the phase offset ϕ is approximately the phase measured at the center of the field; that is, near the origin at

*x*= 0 (or

*t*= 0) (see Eq. 65). Because the central grid field is symmetric about the origin and the phase decreases linearly with position, the phase at grid-field entrance ϕ

_{entry}is related to ϕ (for any spacing

*s*) as follows: With these requirements for ω

_{n}*and ϕ in Equations 68 and 70, respectively, we can write*

_{n}*M*in Equation 65 as follows: Extending Equation 64, the final expression for the modulated place field

_{n}*P*(

_{M}*x*) is as follows: Note that we have used the same synaptic weight function

*A*(

*s*, σ) as for the nonmodulated case in Equation 64. To obtain the membrane potential of the output cell, we perform the convolution between

_{n}*P*(

_{M}*x*) and an EPSP kernel

*D*, top (normalized). The peaks of

*V*

_{out}(

*x*), which are shown in Figure 6

*D*, bottom, exhibit phase precession.

To derive an approximation for the phase-precession range of the output place field, we focus on the fast oscillatory component, which is determined by *M _{n}* in Eqs. 71 and 72. The (angular) spatial frequency of the modulation (coefficient of

*x*in Eq. 71) is as follows: where the term Ω/0.7

*s*denotes the spatial slope of phase recession of a grid cell with spacing

_{n}*s*. Note that because phase decreases with position, the true slope is negative; that is, −Ω/0.7

_{n}*s*.

_{n}The slope of the output place field can be obtained from the mean modulation 〈*M _{n}*〉 with respect to the synaptic weight distribution

*A*(

*s*, σ). We can also estimate the phase-precession slope of the output place field by simply considering the mean spatial slope 〈

_{n}*s*〉 at the input. This approximation corresponds to the first term of a Taylor expansion of 〈

_{n}*M*〉. The phase-precession slope of the output place field is thus Ω/(0.7 · 〈

_{n}*s*〉), where the mean spacing is as follows: The mean slope of phase precession at the output is then Ω/(0.7 · 6.5σ). With the size 3σ of the output place field and the corresponding slope Ω/(0.7 · 6.5σ), we can calculate the phase-precession range Δϕ at the output as follows: Note that the output range Δϕ ≈ 0.66Ω is smaller than the input range Ω, which means that range is reduced in the grid-to-place transformation, and that Δϕ is independent of the size σ of the output place field, which means that slope-size matching is conserved from the MEC grid fields to the hippocampal place fields. For an input range of Ω = 250° (see Fig. 6

_{n}*A–C*), the value of the phase-precession range at the output is Δϕ ≈ 0.66 · 250° ≈ 165°, and this value matches the numerically obtained phase-position plot in Figure 6

*D*, bottom.

## Results

Using a minimal computational model, we investigate inheritance of phase precession; that is, the unidirectional transmission of phase precession from one region to another, for example, in the hippocampal formation. To this end, we consider a population of phase-precessing neurons in one region projecting to a single output cell in another region.

### Phase precession in a CA3-CA1 network model

As a particular example, we first focus on the CA3-CA1 network of the rodent hippocampus. Can phase precession in a CA1 cell be explained by phase precession observed in a population of CA3 cells? To study this case, we consider *N* CA3 place cells with identical place fields that project to a single CA1 cell and model the average firing rate of the CA3 pyramidal cells during the traversal of a place field (Fig. 1*A*). This population firing rate, here as a function of time, reflects the activity of cells that phase precess because this rate oscillates at a frequency (here 8.5 Hz) that is slightly larger than the theta frequency (here 8 Hz). Peak firing rates are about 15 spikes/s in the center of the simulated place field and rates decay to zero outside the field. Using such a firing-rate profile, we generate spike times by means of an inhomogeneous Poisson process. The spike times of *N* = 200 statistically independent, but otherwise identical, simulated CA3 pyramidal cells are shown in the raster plot in Figure 1*B*. Furthermore, each input spike is assumed to elicit an EPSP (Fig. 1*C*, inset) in a target CA1 cell. The depolarization received by a CA1 cell is thus the sum of many EPSPs reflecting a transient excitatory CA3 input (Fig. 1*C*). This contribution to the membrane potential of the simulated CA1 cell is noisy (due to the Poisson process) and is delayed and scaled (due to filtering by the EPSP) with respect to the firing rate of the CA3 cell shown in Figure 1*A*. To simulate the full subthreshold membrane potential of a CA1 cell that receives transient excitatory phase-precessing input from CA3 place cells, we add an ongoing intracellular oscillation at theta frequency, and this oscillation rides on top of a resting membrane potential at −70 mV (Fig. 1*D*). This ongoing intracellular oscillation is assumed to be: (1) typically smaller in amplitude than the maximum amplitude of the excitatory component, (2) present independent of any place-field activity, and (3) phase locked to the extracellularly recorded LFP in the theta band. This ongoing intracellular oscillation could be generated by inhibitory input.

By summing the ongoing inhibitory and the transient excitatory input, we obtain the oscillatory membrane potential of the CA1 model neuron as shown in Figure 1*E*. To determine whether this subthreshold membrane potential reflects phase precession, we investigated the local maxima and assigned an LFP phase to each peak. The peak phases are plotted as a function of time in Figure 1*F*. Within the place field, there is clear phase precession; that is, there is a strong correlation between phase and time, and the slope is negative because the peak phases are mainly determined by the phase-precessing excitatory input from CA3. Outside of the place field, the phase of the peaks is constant at ∼360° because peaks are determined by the inhibitory input that is locked to the LFP theta oscillation.

### Subthreshold signatures of phase precession in CA1 pyramidal cells can be explained by inheritance

The simulated subthreshold membrane potential of a CA1 pyramidal cell (Fig. 1*E*,*F*) exhibits three key features, as observed in *in vivo* whole-cell recordings in awake animals that cross a place field. First, the mean membrane potential rises and falls in a ramp-like manner, which is reminiscent of the CA1 (subthreshold) place field (Lee et al., 2006; Harvey et al., 2009; Epsztein et al., 2011). In our model, the subthreshold place field in CA1 is a result of the CA3 place-field mean activity, which depolarizes the membrane potential. The place field in CA1 is thus inherited from a population of CA3 cells that have overlapping place fields. Second, the amplitude and the frequency of the oscillations of the membrane potential are larger within the place field than outside (Harvey et al., 2009). To explain the larger amplitude using our model, we assumed that the amplitude of the oscillatory contribution of the excitatory input from CA3 is larger than the amplitude of the ongoing inhibitory oscillation, a constraint that we will later use to estimate parameters of the model. To explain the higher frequency, we note that the excitatory contribution originates from CA3 phase-precessing cells that are oscillating at a frequency higher than LFP theta oscillations, and are thus higher than the ongoing inhibitory oscillation, which was assumed to have the same frequency as the theta LFP. Third, peaks of the subthreshold membrane potential show phase precession (Harvey et al., 2009). In the computational model, phase precession in the CA1 membrane potential is inherited from phase precession of the inputs from CA3. Overall, the simulation results suggest that CA3 excitatory input and an ongoing inhibitory oscillation are sufficient to explain the main subthreshold features of the membrane potential of CA1 pyramidal cells that show phase precession (Harvey et al., 2009).

### Parametrization of the computational model

How are the subthreshold features of phase precession related to the parameters of the computational model of a CA3-CA1 network? Let us first briefly describe the model in mathematical terms (for details, see Materials and Methods). The input from CA3 is described by the time-dependent firing rate λ(*t*) = λ_{0}[1 + *C* cos(2π*f*_{λ}*t* − ϕ_{λ})]*x*(*t*) (see also Eq. 4) with mean firing rate λ_{0} in the center of the place field, modulation depth *C* of the oscillatory component of the firing-rate oscillation at frequency *f*_{λ} and phase ϕ_{λ}. The Gaussian function *x*(*t*), with maximum value 1 and size σ, describes the temporal extent of a CA3 place field (Fig. 1*A*). This firing-rate model is used to generate spikes, each spike evokes an EPSP, and EPSPs are linearly integrated. An EPSP is described by an α function with amplitude ϵ_{max} and time constant τ (Fig. 1*C*, inset). A further important model parameter is the number *N* of CA3 cells that provide input to a single CA1 cell. The ongoing theta oscillation of the CA1 membrane voltage (Fig. 1*D*) is quantified by the amplitude *B*, the theta frequency *f*_{θ}, and the phase ϕ_{θ} (with respect to the extracellular theta LFP): *V*_{θ}(*t*) = *B*[cos(2π*f*_{θ}*t* − ϕ_{θ}) − 1] (see also Eq. 1).

Some model parameters are constrained by experimental data (mean ± SE): λ_{0} = 12.4 ± 4 spikes/s (estimated from main text and Fig. 7*H* in Mizuseki et al., 2012), *f*_{λ} = 8.6 ± 0.3 Hz (estimated from Fig. 14*A* in Mizuseki et al., 2012), ϕ_{λ} = 190 ± 30° (estimated from Fig. 2*D* in Harris et al., 2002), *f*_{θ} = 8.1 ± 0.2 Hz (estimated from Figs. 1⇓⇓⇓–5 in Geisler et al., 2010), *B* = 0.7 ± 0.1 mV (estimated from Fig. 5*C* in Harvey et al., 2009), and τ = 10 ± 3 ms (estimated from Fig. 3*D* in Magee and Cook, 2000). Note that τ ≈ τ* _{d}*/2 where τ

*is the decay constant in Magee and Cook (2000). However, other parameters are not easily accessible. In particular, the phase ϕ*

_{d}_{θ}of ongoing intracellular theta oscillations outside of the place field (with respect to theta LFP recorded in the CA1 pyramidal layer) is unclear in awake, behaving animals. Furthermore, for the input from CA3 to CA1, the modulation depth

*C*of the population firing rate is not available. Finally, the number

*N*of CA3 cells that drive one CA1 cell and the

*in vivo*EPSP amplitude ϵ

_{max}are, so far, free parameters of the model. In what follows, we constrain these parameters by comparing the computational model with available data on the properties of place fields and phase precession, including an analysis of the noise in the membrane potential.

### Subthreshold features of phase precession constrain model parameters

While an animal traverses a place field, CA3 neurons fire sequences of spikes, and the spikes evoke EPSPs that contribute to the membrane voltage of a CA1 neuron. This sum of EPSPs is variable in each run because the activity of CA3 neurons is variable. To account for such variability in the inheritance model, spikes are generated in each CA3 neuron by means of an inhomogeneous Poisson process. The variability of the membrane voltage due to such “shot noise” critically depends on the number *N* of CA3 cells that generate the place-field response in the CA1 cell. Intuitively, we expect a more faithful representation of the CA3 population activity in a CA1 cell for an increasing number *N* of CA3 neurons. In Figure 2, we show results of numerical simulations of CA1 membrane potential traces for values of *N* between 30 and 260. Note that, here, we only consider the phase precessing—that is, excitatory—contribution to the membrane potential. To allow for an easy comparison of membrane-potential traces associated with values of *N* that span one order of magnitude, we scale EPSP amplitudes proportional to 1/*N*. The appearance of a simulated CA1 membrane-voltage trace also depends on the modulation depth *C* of the oscillatory firing rate of the simulated CA3 cells. Figure 2 shows that increasing *C* from 0.3 to 0.7 (for fixed *N*) enhances the oscillatory component of voltage traces.

The simulated voltage traces in Figure 2 show that the model parameters *N* and *C* have a major impact on the shape of the CA1 membrane potential. However, there is a trade-off between *N* and *C*. Indeed, the oscillations of voltage traces on the diagonal from top-left to bottom-right in Figure 2 are similar. To quantify differences and similarities in the appearance of voltage traces, which is essential to constrain the ranges of the model parameters *N* and *C*, we introduce a signal-to-noise ratio. The signal-to-noise ratio ρ is defined as the average amplitude of the membrane potential oscillation divided by the average amplitude of the shot noise in the center of the place field (for details, see Materials and Methods, Eq. 29):
We note that ρ is independent of the EPSP amplitude, which justifies the voltage scaling of the traces in Figure 2. The larger the ρ, the greater the amplitude of the oscillations compared with the shot noise or, more qualitatively, the clearer the oscillatory structure of the CA1 voltage trace. For example, the voltage trace for *N* = 30 and *C* = 0.3 shows a membrane-potential trace for which the oscillations are barely distinguishable from the noise, which is described by a rather low value ρ = 0.7. In Figure 2, the signal-to-noise ratio ρ increases from bottom left (ρ = 0.7) to top right (ρ = 6.5) because both *N* and *C* are increased. On the diagonal from top left to bottom right, the signal-to-noise ratio is constant (ρ = 2.2) because of the particular choices of *N* and *C*. These traces are most similar to the CA1 voltage traces recorded by Harvey et al. (2009) in awake, behaving mice.

So far, the signal-to-noise ratio suggested by experimental data (e.g., ρ = 2.2) constrains the model parameters *N* and *C*, but does not lead to particular values. Moreover, the EPSP amplitude ϵ_{max} cannot yet be determined. We therefore need to relate further properties of experimental subthreshold membrane voltage traces to the computational model. Characteristic quantities that can be extracted from the data are, for example, the oscillation amplitude Δ*V*_{osc} and the mean depolarization ramp Δ*V*_{ramp} in the center of the field. Harvey et al. (2009) found that Δ*V*_{osc} = 1.3 ± 0.4 mV and Δ*V*_{ramp} = 2.7 ± 0.4 mV.

A thorough theoretical analysis of our inheritance model reveals a direct relationship between model parameters and the oscillation amplitude Δ*V*_{osc} and the mean depolarization ramp Δ*V*_{ramp} (for details, see Materials and Methods, Eqs. 16 and 17):
and
where e = exp(1) ≈ 2.71828 is Euler's number. Both equations are excellent approximations if the amplitude *B* of ongoing oscillations is smaller than the oscillation amplitude Δ*V*_{osc} of the transient component in the center of the field.

Equation 77 for ρ, Equation 78 for Δ*V*_{osc}, and Equation 79 for Δ*V*_{ramp} are sufficient to determine the three free model parameters *C*, *N*, and ϵ_{max} (for details, see Materials and Methods, Eqs. 31, 32, and 33):
Furthermore, we use standard rules of error propagation applied to Equations 80, 81, and 82 and the experimental values for τ, Δ*V*_{osc}, Δ*V*_{ramp}, *f*_{λ}, λ_{0}, and ρ to estimate the mean and error (i.e., SEM) of model parameters *C*, *N*, and ϵ_{max} as follows:
To summarize, using the computational model and salient features of subthreshold voltage data, we are able to explicitly predict values for the model parameters *C*, *N*, and ϵ_{max}. The value for ϵ_{max} is comparable to the values of EPSPs (≈0.2 mV) recorded *in vitro* (Magee and Cook, 2000). An estimate for the modulation depth *C* from the literature requires recordings of the activity of a population of CA3 cells with overlapping place fields during a single place-field traversal. Such data are not available. A preliminary estimate, however, can be obtained from the averaged (across multiple traversals) firing rates of CA1 cells. From Skaggs et al. (1996, their Fig. 3*A*), we estimated *C* ≈ 0.5, which matches our prediction in Equation 83. We note that the prediction in Equation 84 for the number *N* of presynaptic CA3 cells constitutes a lower bound. To arrive at this estimate, we neglected other sources of noise and variability and also correlations in the input, which typically increase *N* (see Discussion for details).

### Phase precession in a bursting model neuron

We have shown that the inheritance model can account for subthreshold features of phase precession. Phase precession, however, is also a suprathreshold phenomenon in which spikes of place cells shift to earlier phases relative to theta LFP. Therefore, in the context of a spiking model neuron, we investigated whether the input from CA3 cells that exhibit phase precession is sufficient to produce phase precession of spikes of a CA1 cell. For this purpose, we use a two-compartment model to: (1) segregate inhibitory input onto the perisomatic compartment from excitatory input onto the dendritic compartment, (2) include effects of active conductances in dendrites, and (3) account for bursting. We therefore use the model developed by Pinsky and Rinzel (1994) for CA3 neurons and adapted for CA1 neurons by Kamondi et al. (1998). A schematic is shown in Fig. 3*A* (for details, see Materials and Methods).

The two inputs to the CA1 model neuron are dendritic (*I _{d}*) and somatic (

*I*) current densities (μA/cm

_{s}^{2}), referred to hereafter simply as currents. The dendritic current

*I*represents the total phase-precessing signal coming from CA3. More precisely, each input spike elicits an EPSC. The current

_{d}*I*is thus a sum of EPSCs and this summed current is applied to the dendritic compartment of the model neuron. The somatic current

_{d}*I*represents an inhibitory oscillation that is phase locked to the LFP (Kamondi et al., 1998) and is applied to the somatic compartment.

_{s}Figure 3*B* shows an example of the time courses of both currents and the CA1 membrane potential for one traversal of the place field. The two-compartment model generates isolated action potentials and bursts, and the spiking activity shows phase precession, here with respect to time, similar to CA1 pyramidal cells *in vivo* (O'Keefe and Recce, 1993; Pastalkova et al., 2008; Schmidt et al., 2009). Simulated membrane potentials for several traversals through the CA1 virtual place field are shown in Figure 3*C* and the pooled spike phase versus time plot in Figure 3*D* again demonstrates phase precession.

In the spiking model, the intrinsic ionic currents in the dendrite and soma are responsible for the bursting dynamics. However, the phase precession in the CA1 model cell is essentially due to the input from CA3 phase-precessing cells. Indeed, the spikes and bursts occur at the maxima of the underlying oscillations of the membrane potential, as shown in Figure 3*C*. This means that the subthreshold membrane potential peaks are a good predictor of suprathreshold behavior, which is consistent with data from Harvey et al. (2009), Domnisoru et al. (2013), and Schmidt-Hieber and Häusser (2013). Because of this tight link between subthreshold and suprathreshold phase precession in our model, hereafter we restrict the analyses to the subthreshold dynamics only.

### Intracellular theta oscillations modulate phase precession of CA1 model cells

In the previous sections, we examined in detail the role of excitatory input from CA3 phase-precessing cells but largely neglected the influence of the ongoing intracellular theta oscillations in CA1, which are phase locked to the extracellular theta LFP. Now, we elucidate the role of these ongoing oscillations, which may be mediated by somatic inhibitory input into CA1 pyramids. To this end, we vary the phase ϕ_{θ} and the amplitude *B* of the ongoing theta oscillations in our model while the excitatory input is kept fixed (Fig. 4). We restrict the analysis to a mean-field model; that is, we neglect the shot noise. The excitatory input is therefore the firing rate of a population of CA3 cells convolved with an EPSP.

To evaluate the influence of an ongoing, putatively inhibitory, theta oscillation of the membrane potential on phase precession in CA1, we compare this more involved case with the simpler case without such ongoing oscillations (Fig. 4, red); we note that there is clear phase precession in CA1 in the case without ongoing oscillations.

We begin by studying the impact of the phase ϕ_{θ} on phase precession while the amplitude *B* = 1 mV is fixed (Fig. 4,*A*,*B*). To demonstrate how the range and the overall shape of phase precession depends on the phase ϕ_{θ} of ongoing somatic inhibitory input, we first consider the case of ϕ_{θ} ≈ 0 = 360°. In this case, the peak phases of the ongoing somatic oscillation are in phase with the peaks of the extracellular LFP. Interestingly, for ϕ_{θ} ≈ 0°, the slope, range, and entry phase of phase precession are larger than for the standard case with no ongoing theta oscillations. For ϕ_{θ} ≈ 0°, the phases of the CA1 membrane voltage peaks at field entry are ≈300° (Skaggs et al., 1996; Mizuseki et al., 2009, 2012), which can be explained by the entrance phase of CA3 firing (≈200°, Harris et al., 2002), a subsequent phase delay due to synaptic filtering of the EPSP (≈50° = 1.8τ · *f*_{θ} · 360°; Eq. 14), and a delay due to the inhibitory oscillation (≈50°, estimated from Fig. 4*B*). Oscillatory inhibitory input therefore can explain why the range of phase precession in CA1 is larger than in CA3 (Harris et al., 2002; Mizuseki et al., 2009, 2012). Overall, the case ϕ_{θ} ≈ 0° matches *in vivo* data well. In contrast, for other phases of the ongoing somatic oscillation, for example, ϕ_{θ} ≈ 120° and ϕ_{θ} ≈ 240°, phase precession does not match *in vivo* data: the slope and the range of phase precession are reduced and, at field entry, the phase is smaller compared with the standard case of pure CA3 excitatory input.

To study the influence of the amplitude of the ongoing intracellular theta oscillations on phase precession, we fix the phase at the value ϕ_{θ} ≈ 0°, which matches *in vivo* data on phase precession best, and vary the amplitude *B* (Fig. 4*C*,*D*). Phase precession strongly depends on the amplitude *B* of ongoing somatic inhibitory input. The range of phase precession is largest for moderate values of *B*, for example, *B* = 1 mV (as in Fig. 4*A*,*B*) and *B* = 2 mV. In contrast, for *B* = 5 mV, phase precession practically disappears because the ongoing intracellular oscillation at theta frequency dominates the transient excitatory component. Conversely, for *B* = 0, that is, excitation only, and 0.5 mV, the range of phase precession is reduced compared with the case *B* = 1 mV.

The detailed phase relation of the traces for *B* = 0 and *B* = 1 mV in Figure 4*D* is even consistent with data from Royer et al. (2012): for ϕ_{θ} ≈ 0°, a reduction of the amplitude of inhibition leads to smaller phases near field entry (phase advance), whereas toward field exit phases are larger (phase delay). The assertion that ϕ_{θ} ≈ 0° is also consistent with the following data. In awake animals during the theta state, parvalbumin-expressing (PV) interneurons preferentially fire before the trough of the theta LFP, at ∼130° (Royer et al., 2012; Varga et al., 2012), and project to the perisomatic region of CA1 pyramidal neurons. Because of the inhibitory nature of this projection and a small delay mainly due to synaptic filtering, the evoked oscillation of the membrane potential in CA1 pyramids can be assumed to have minima at ∼180° and, therefore, maxima at ∼0°, as predicted.

In summary, the inheritance model of phase precession predicts that the ongoing intracellular theta oscillation in CA1 pyramidal cells is not necessary for phase precession. Phase precession remains without such intracellular theta oscillations because phase precession is inherited from the excitatory input from CA3. However, the interaction of a specific ongoing oscillation, possibly due to inhibitory interneurons, and transient excitatory oscillations from CA3 pyramids can modulate phase precession in CA1 in a characteristic way. For ϕ_{θ} ≈ 0°, the range and entry phase of phase precession are reduced if inhibitory input is suppressed, which is supported by data from Royer et al. (2012). Therefore, the inheritance model of phase precession predicts that the ongoing intracellular theta oscillation in CA1 pyramidal cells is in phase with the extracellular LFP theta oscillation; that is, ϕ_{θ} ≈ 0°.

### Inheritance explains phase precession for a variety of place-selective responses

So far, we have examined a CA1 place cell that inherits phase precession from a subset of CA3 cells that share a common place field; that is, the input corresponded to identical place fields of many different CA3 cells. We now relax this strong assumption and study a more general case of inheritance in which the place fields at the input are allowed to be spatially distributed. This case is equivalent to place fields at the input that are distributed uniformly but the respective synaptic weights in the target cell are distributed spatially (Malhotra et al., 2012).

To study spatially distributed inputs in the inheritance model, we still assume that each input cell shows phase precession within its place field; that is, input place fields are characterized by a common frequency *f*_{λ} of modulatory activity (*f*_{λ} > *f*_{θ}). Moreover, input place fields all have the same size, which is parameterized by the width σ. Furthermore, pairs of overlapping input place fields are aligned such that the theta-time scale delay between firing-rate peaks is correlated to the distance between place-field centers (Dragoi and Buzsáki, 2006; Geisler et al., 2010; for details see Materials and Methods). Here, we neglect the ongoing inhibitory oscillations, which were argued to be modulatory. As in earlier sections, we assume that an output neuron linearly integrates the excitatory input, producing a somatic voltage as a response.

To connect the analysis of phase precession for spatially distributed place-field centers to analyses already discussed, we first consider again the case of identical input place-field centers (Fig. 5*A*, top); that is, all input place fields are perfectly overlapping. This simplified scenario corresponds to the topology of our previous results regarding the CA3-CA1 system and serves as a reference for the other distributions of input fields in Figure 5. By summing the depolarization provided by the identical input place fields, we obtain an output place field (Fig. 5*A*, middle) that is essentially the same as the input place field. As expected, there is phase precession throughout the output place field, as evidenced by the negative slope of the corresponding phase-position plot (Fig. 5*A*, bottom).

Next, we consider a Gaussian distribution of the centers of input place fields; the width of this distribution is characterized by σ* _{d}* (Fig. 5

*B*). Adding such distributed input place fields, we obtain an output place field that is larger than each individual input place field (as in Fig. 5

*A*), and the size

*of the distribution of field centers (Eq. 51). This scenario can explain why, on average, the size of place fields increases from CA3 to CA1 and to the subiculum (Kim et al., 2012; Mizuseki et al., 2012). Moreover, the slope of the phase-position plot at the output is shallower than at the input, although the range of phase precession is similar. Therefore, the inheritance model can explain why the phase-precession slope is inversely proportional to the field size (Huxter et al., 2003; Terrazas et al., 2005; Dragoi and Buzsáki, 2006; Kim et al., 2012).*

_{d}Another interesting case is a spatially uniform distribution of the place-field centers (Fig. 5*C*). This case might correspond, for example, to the input of a basket cell in CA1. Adding such uniformly distributed inputs, we observe an oscillation of the membrane potential of the output cell, but there is no spatial preference. The corresponding phase-position plot indicates that peaks of the membrane potential do not show phase precession. The phase of the peaks of the membrane potential oscillation is at ∼130° and the membrane potential is oscillating precisely at theta frequency *f*_{θ}, although the phase-precessing input cells were oscillating at a frequency *f*_{λ} > *f*_{θ} (see also Geisler et al., 2010). This behavior of the output cell is reminiscent of the activity of basket cells in the hippocampus, which receive excitatory input from place cells. Basket cells fire phase locked to theta, the preferred firing phase is ∼130° (Varga et al., 2012), and there is little spatial preference (Frank et al., 2001).

Next, we consider a ramp-like distribution of place-field centers (Fig. 5*D*) in which the density of place fields at the input (or, equivalently, the synaptic weights) increases linearly along a specified axis. The resulting output place field is, not surprisingly, also ramp like: the mean depolarization increases steadily and falls off abruptly. Furthermore, this ramp-like depolarization is modulated by an oscillation that shows a characteristic shape of phase precession, with an initially shallow slope, some plateau <180°, and, finally, some steeper phase precession. Such ramp-like neurons are found in the ventral striatum (van der Meer and Redish, 2011) and, remarkably, the observed shape of the phase precession resembles the shape generated by our model.

Finally, we test whether hippocampal phase precession can be inherited from the MEC. Principal cells in the MEC fire in a grid-like manner in open environments (Hafting et al., 2005) and exhibit phase precession, even after inactivation of the hippocampus (Hafting et al., 2008). Furthermore, grid cells project to the hippocampus (Zhang et al., 2013). Is it possible to obtain a phase-precessing place cell from a population of phase-precessing grid cells? Analogous to our derivations concerning place-field distributions, we assume an input layer of phase-precessing grid cells projecting to a single output hippocampal cell. In the context of our previous results regarding inheritance of phase precession to CA1, the input grid cells could correspond to cells from MEC layer III that project to a cell in CA1.

In the one-dimensional case we study here, the grid-field activity is periodic along the axis where the virtual rat is running (Brun et al., 2008). Following our previous approach for CA3-CA1 inheritance (Fig. 1), we implemented phase precession within grid fields through a time-dependent firing-rate modulation (Fig. 6*A–C*). In our simulations, grid cells have different field spacings, and the size of the firing field, which is defined by a 20% threshold (Brun et al., 2008), is ∼0.7 times the grid spacing. Grid cells also have phase-precession slopes that are inversely proportional to the field size, with a constant phase-precession range of ∼250° (Hafting et al., 2008; Reifenstein et al., 2012a) and the spike phase at field entry is ∼200° (Climer et al., 2013; Reifenstein et al., 2012b).

The synaptic weights of the projection from MEC layer III to CA1 are calculated as Fourier coefficients such that a linear combination of periodic grid fields produces a nonperiodic place field at the output (for details, see Materials and Methods and Solstad et al., 2006 for the original derivation in two dimensions). The resulting synaptic weight distribution in Figure 6*E* is asymmetric with a steeper rise and shallower decay. For the inheritance model, we used a uniform distribution of grid spacings, although recent studies (Stensola et al., 2012; Kitamura et al., 2014; Ray et al., 2014) have highlighted the existence of discrete modules or clusters in the MEC that reflect a degree of anatomical spatial organization. Including the effect of such clusters in our model would result in a coarser discretization of the grid spacing of the synaptic weight distribution without affecting the main results.

The membrane potential of the output cell reflects a single place field of width 3σ = 0.66 m (Fig. 6*D*), a value that is close to the smallest grid spacing *s* with a significant synaptic weight (Fig. 6*E*). The membrane potential peaks of the output place field exhibit phase precession (Fig. 6*D*, bottom); however, the range of phase precession is only ∼145°. This small range, which is the product of the size and the phase-precession slope of the output field, can be explained by the grid-to-place inheritance. The size ∼3σ of the output field is related to the shape, for example, the steep rise, of the synaptic weight distribution in Figure 6*E* (solid arrow). The phase-precession slope of the output field, however, is largely determined by the mean slope of the input grid fields and this mean slope is related to mean grid spacing of the input population. Figure 6*E* (dashed arrow) indicates that the mean grid spacing of the input population is ∼6.5σ = 1.4 m (Eq. 75) and therefore the mean grid-field size is ∼0.7 · 1.4 m = 1 m. Accordingly, the mean slope of phase precession in the input is ∼250°/1 m = 250°/m. Using this slope and the obtained size of the output place field (3σ = 0.66 m), we estimated the range of phase precession as 250°/m · 0.66 m = 165° (see also Materials and Methods, Eq. 76). This estimation assumed a linear phase-position relation, although the phase-position relation in Figure 6*D* is actually sigmoidal. Therefore, the numerically obtained phase range of ∼145° is smaller. The saturation of the phase at the borders of the field indicates the transition to a constant phase of membrane-voltage peaks outside the place field, a behavior that is similar to Figure 5*C*.

The phase of the first spike (entry phase) in the modeled output place field in Figure 6*D* is ∼200°, which is equal to the entry phase 200° of the input grid fields. This seemingly paradox zero phase difference, however, must not be interpreted as a zero transmission delay between the two regions, but rather can be explained by the cancellation of two effects. The ∼100° reduction of the phase range from input to output leads to an ∼50° advance of the onset phase because of the symmetry of the place field around its center. Conversely, there is a phase delay due to synaptic filtering of the EPSP and axonal transmission, which is also ∼50° (≈1.8τ · *f*_{θ} · 360°; see Eq. 14).

Phase precession can thus be propagated from MEC layer III to CA1. Interestingly, the resulting phase-position relations due to input from MEC or from CA3 significantly overlap within a theta cycle (cf. Fig. 5*A* for CA3, Fig. 6*D* for MEC layer III). Therefore, our model predicts that both MEC layer III and CA3 can contribute to the phase precession observed in CA1. More generally, a grid-to-place transformation is sufficient to account for inheritance of phase precession from entorhinal cortex to the hippocampus and can explain the observed reduction of the range of phase precession (Harris et al., 2002; Mizuseki et al., 2009; Schmidt et al., 2009; Reifenstein et al., 2012a).

Overall, inheritance of phase precession can explain a variety of subthreshold and suprathreshold responses in the hippocampus and related structures in terms of place selectivity and presence of phase precession.

## Discussion

We developed a computational model that explains inheritance of phase precession. Our main findings are as follows: (1) phase precession in a small subset of CA3 or MEC principal cells is sufficient to explain phase precession in a CA1 pyramidal cell; (2) inhibitory input to the CA1 cell is not necessary to generate phase precession, but inhibition can increase the slope and enlarge the range of phase precession; and (3) the spatial distribution of field centers and the width of fields of the input population determine not only the spatial selectivity, but also the phase-precession characteristics of the target cell.

According to this inheritance model, the subthreshold signatures of phase precession in a CA1 cell (Harvey et al., 2009) reflect the firing-rate characteristics of phase-precessing input cells. Particularly, the mean depolarization (≈4 mV) of a CA1 cell in the center of a place field (i.e., the subthreshold place field) is the result of increased excitatory drive; the oscillatory component of the CA1 membrane potential (≈2 mV amplitude in the field's center) is due to phase-precessing feedforward input; and voltage traces are noisy due to variable spiking input, but the oscillation amplitude in the theta frequency range is larger than the SD of the noise in single runs through a place field. Finally, we showed that peaks of the modeled subthreshold oscillations show phase precession (Harvey et al., 2009; Domnisoru et al., 2013; Schmidt-Hieber and Häusser, 2013) and, using a spiking model neuron, we confirmed that the peaks predict phase precession of action potentials.

To generate phase precession in a CA1 pyramidal cell, the experimental constraints outlined above imply that we need *N* > 200 active input cells, with a population firing rate with a modulation depth *C* ≳ 0.6. To arrive at these predictions, we assumed that peak firing rates of CA3 cells are the same although they are variable (Leutgeb et al., 2006) and we assumed identical synaptic weights and transmission probabilities for individual projections. Furthermore, our CA3 population firing-rate model did not include correlations that arise, for example, from task-dependent switching of reference frames (Jackson and Redish, 2007) or from assembly-level organization (Foster and Wilson, 2007). Such correlations hint at the importance of network-level interactions during phase precession (Lisman and Redish, 2009). Such additional variability and correlations of the input, which are not included in our Poisson model, typically increase the noise in the membrane voltage of a CA1 cell. Our model thus underestimates the true values of *N* and *C*, and our predictions are therefore lower bounds. Importantly, our CA3 population firing-rate model is not specific to any theory of how phase precession is generated *de novo*.

We simulated the inhibitory input to a CA1 cell by an ongoing oscillation of the membrane potential that is coherent with the theta LFP but is independent of place-specific activity (Harvey et al., 2009; Epsztein et al., 2011). In our model, this oscillation was parametrized by the phase ϕ_{θ} of peaks with respect to the theta-LFP peaks and by its amplitude *B*. Only for ϕ_{θ} ≈ 0 and amplitudes *B* ≲ 1 mV, which are smaller but comparable in magnitude to the maximum oscillation amplitude (≈2 mV) due to the incoming excitation, does the inhibition increase the range and the slope of phase precession of a CA1 cell compared with the case in which the CA1 cell receives excitation only. That inhibition only modulates phase precession, but is not necessary for its generation, receives support from Royer et al. (2012), who showed that transiently silencing PV interneurons reduced range and slope of phase precession of neighboring pyramidal cells in a characteristic way, which is reproduced in detail by the inheritance model.

The specific prediction ϕ_{θ} ≈ 0 of the inheritance model distinguishes it, for example, from the somatodendritic interference model (SDI) of phase precession (Kamondi et al., 1998; Magee, 2001; Harris et al., 2002; Mehta et al., 2002), which requires maxima of the inhibition at ϕ_{θ} ≈ 180°. In the SDI model, phase precession within a single neuron arises from the interference between a slow, ramp-like excitatory input to dendrites and an oscillating, inhibitory input to the soma. Phases of spikes of a CA1 pyramidal cell at the entrance of a place field then occur at ≈ 0° because, at this phase, the excitation just exceeds the minima of the inhibitory oscillation. Spike phases then precess because of increasing excitatory input. For decreasing excitatory input toward the end of the place field, some additional adaptation of CA1 cells or asymmetric excitation suppresses phase “recession.” The SDI model is therefore in stark contrast to the inheritance model. Experiments that determine the theta phase of oscillations evoked by inhibitory input *in vivo* could reject one of them.

In general, a difference between the inheritance model and previous models of phase precession is that the latter can explain how phase precession is generated ab initio. Network models, for example, require considerable recurrent connectivity within the local network (Jensen and Lisman, 1996; Tsodyks et al., 1996; Wallenstein and Hasselmo, 1997; Navratilova et al., 2012). These models are therefore more consistent with the denser recurrent connectivity of CA3 or EC than that of CA1 (Amaral and Witter, 1989; Couey et al., 2013; Pastoll et al., 2013). However, even in CA3 recurrent connectivity may not be necessary to explain phase precession, which could alternatively arise from facilitation of the mossy fiber synapse (Thurley et al., 2008) or inheritance from the dentate gyrus and/or EC (Molter and Yamaguchi, 2008).

We showed that CA3 input is sufficient to generate phase precession in CA1, but CA1 receives excitatory input also from layer III of the MEC, which can generate place-specific activity in CA1 (Brun et al., 2002; Nakashiba et al., 2008). Furthermore, cells in MEC layer III exhibit phase precession (Hafting et al., 2008; Reifenstein et al., 2012b; Climer et al., 2013; but see Mizuseki et al., 2009). We showed that the contributions from CA3 and from MEC layer III generate similar phase-position relations, so both pathways can contribute to the phase precession observed in CA1. Once the mechanism(s) for phase precession in CA3 and MEC are established, they will be sufficient to explain phase precession in CA1 through inheritance (Skaggs et al., 1996; Yamaguchi et al., 1998, 2007).

The specific results of the CA3-CA1 inheritance model can be generalized to other pathways. We found that the spatial distribution of centers of phase-precessing cells of an input population determines not only the place-field selectivity, but also the phase-precession characteristics of an output neuron. Therefore, the inheritance model can account for several further features of phase precession (Figs. 5, 6). First, the phase-precession slope matches the size of the place field so that the range remains approximately constant (Huxter et al., 2003; Dragoi and Buzsáki, 2006). Second, place-field sizes in the input population are smaller than those downstream (Jung and McNaughton, 1993; Kim et al., 2012; Mizuseki et al., 2012; but see Lee et al., 2004). Third, putative basket cells show little spatial preference and have a preferred firing phase of ∼130° (Skaggs et al., 1996; Lapray et al., 2012; Royer et al., 2012; Varga et al., 2012), which can be explained by inputs from place cells with place-field centers that are evenly distributed in space (Geisler et al., 2010). Some weak phase precession in interneurons (Maurer et al., 2006a; Ego-Stengel and Wilson, 2007; Geisler et al., 2007) may be due to spatially uneven excitatory input. Fourth, phase precession in the ventral striatum (van der Meer and Redish, 2011) has a characteristic shape and there is a ramp-like activity of neurons, which is explained by inheritance from hippocampal place fields with increasing density of fields when the animal is approaching a goal location (Hetherington and Shapiro, 1997). Fifth, CA1 cells can have two or more overlapping place fields and, in each subfield, there is phase precession (Maurer et al., 2006b). This spatial overlap of phase precession in a single CA1 cell can be explained by two or more sufficiently synchronously phase-precessing input assemblies (see also Jones and Wilson, 2005, for phase precession in the prefrontal-cortex). Finally, the “grid-to-place” transformation (McNaughton et al., 2006; Rolls et al., 2006; Solstad et al., 2006; Gorchetchnikov and Grossberg, 2007; Blair et al., 2008; de Almeida et al., 2009; Savelli and Knierim, 2010; Cheng and Frank, 2011) can also account for the inheritance of phase precession.

In summary, our results predict that the spatial distribution of field centers (or, equivalently, a corresponding distribution of synaptic weights; Malhotra et al., 2012) and the field widths of a population of phase-precessing cells determine the phase precession of a downstream cell that receives this input. Such distributions may be shaped through synaptic plasticity, which is an interesting topic for future theoretical and experimental work.

To conclude, phase precession can be a result of network, cellular, and inheritance effects combined. Inheritance in particular might help to explain why phase precession is observed in many parts of the brain.

## Footnotes

This work was supported by the Bernstein Center for Computational Neuroscience Berlin (Grant 01GQ1001A), the Bernstein Focus “Neuronal Foundations of Learning” (Grant 01GQ0972 to R.K.), the Deutsche Forschungsgemeinschaft (Grant GRK1589/1), and the NeuroCure excellence cluster. We thank Kay Thurley, Eric Reifenstein, Tiziano D'Albis, Dietmar Schmitz, and Bruce McNaughton for valuable discussions and/or comments on the manuscript.

The authors declare no competing financial interests.

- Correspondence should be addressed to Richard Kempter, Department of Biology, Institute for Theoretical Biology, Humboldt-Universität zu Berlin, Invalidenstrasse 43, 10115 Berlin, Germany. r.kempter{at}biologie.hu-berlin.de