Advantages and detection of phase coding in the absence of rhythmicity

Abstract The encoding of information in spike phase relative to local field potential (LFP) oscillations offers several theoretical advantages over equivalent firing rate codes. One notable example is provided by place and grid cells in the rodent hippocampal formation, which exhibit phase precession—firing at progressively earlier phases of the 6–12 Hz movement‐related theta rhythm as their spatial firing fields are traversed. It is often assumed that such phase coding relies on a high amplitude baseline oscillation with relatively constant frequency. However, sustained oscillations with fixed frequency are generally absent in LFP and spike train recordings from the human brain. Hence, we examine phase coding relative to LFP signals with broadband low‐frequency (2–20 Hz) power but without regular rhythmicity. We simulate a population of grid cells that exhibit phase precession against a baseline oscillation recorded from depth electrodes in human hippocampus. We show that this allows grid cell firing patterns to multiplex information about location, running speed and movement direction, alongside an arbitrary fourth variable encoded in LFP frequency. This is of particular importance given recent demonstrations that movement direction, which is essential for path integration, cannot be recovered from head direction cell firing rates. In addition, we investigate how firing phase might reduce errors in decoded location, including those arising from differences in firing rate across grid fields. Finally, we describe analytical methods that can identify phase coding in the absence of high amplitude LFP oscillations with approximately constant frequency, as in single unit recordings from the human brain and consistent with recent data from the flying bat. We note that these methods could also be used to detect phase coding outside of the spatial domain, and that multi‐unit activity can substitute for the LFP signal. In summary, we demonstrate that the computational advantages offered by phase coding are not contingent on, and can be detected without, regular rhythmicity in neural activity.


| INTRODUCTION
In the central nervous system, phase coding refers to the encoding of information in the phase of neuronal activity with respect to an ongoing oscillation in the local field potential (LFP). Theoretically, the coding of information in spike phase offers numerous advantages over a pure rate code. Phase coding can be used to facilitate rapid pattern classification (Thorpe, Delorme, & Van Rullen, 2001) and incurs a lower metabolic cost than an equivalent spike rate code (Fries, Nikolic, & Singer, 2007). It may also be useful for the multiplexing of information (Panzeri, Brunel, Logothetis, & Kayser, 2010); disambiguating stimuli that generate similar spike rates (Kayser, Montemurro, Logothetis, & Panzeri, 2009;Montemurro, Rasch, Murayama, Logothetis, & Panzeri, 2008); and preventing interference between simultaneously presented stimuli (Fries et al., 2007;Jensen et al., 2014).
To date, it has often been assumed that the absence of a baseline oscillation with relatively constant frequency and high amplitude, characterized by a narrow peak in the LFP power spectrum or spike train temporal auto-correlogram, is evidence for a lack of phase coding (Yartsev, Witter, & Ulanovsky, 2011; but see Barry, Bush, O'Keefe, & Burgess, 2012). Nonetheless, firing phase is independent of frequency, such that phase coding can be implemented with a baseline oscillation that varies dynamically over a wide range of frequencies (Blair, Wu, & Cong, 2013;Orchard, 2015). Indeed, it has recently been demonstrated that place cells in bat hippocampus exhibit both phase locking and phase precession relative to aperiodic, low frequency fluctuations in the LFP (Bush & Burgess, 2019;Eliav et al., 2018). Experimental evidence also hints at such a scheme in the human brain, where narrow peaks in power spectra are rare but spiketriggered LFP averages indicate that neural firing is phase locked to ongoing oscillations (Jacobs et al., 2013;Jacobs, Kahana, Ekstrom, & Fried, 2007). This is supported by behavioral evidence showing that increased spike phase coherence correlates with mnemonic performance, even in the absence of an associated peak in the LFP power spectrum (Rutishauser, Ross, Mamelak, & Schuman, 2010).
Here, we demonstrate that spike phase can robustly encode information in the absence of a constant frequency, high amplitude baseline oscillation. As an example, we examine the phase code for location exhibited by grid cells of the medial entorhinal cortex. Using LFP recordings from depth electrodes in the hippocampal formation of pre-surgical epilepsy patients performing a spatial memory task (Bush et al., 2017), we demonstrate that phase precession in simulated grid cells can be robustly maintained in the absence of any clear peak in the LFP power spectra or spike train temporal auto-correlogram, consistent with recent data from flying bats (Eliav et al., 2018). We then demonstrate that this allows simulated grid cells to multiplex information beyond that encoded in firing rate alone (Fiete, Burak, & Brookings, 2008;Mathis, Herz, & Stemmler, 2012), with firing phase indicating movement direction analogous to experimentally observed "theta sequences" of activity along an animal's current trajectory (Burgess, Recce, & O'Keefe, 1994;Johnson & Redish, 2007;Skaggs, McNaughton, Wilson, & Barnes, 1996;Zutshi, Leutgeb, & Leutgeb, 2017). Extending previous observations in place cells (Jensen & Lisman, 2000), we also show that firing phase can be used to improve the accuracy of decoding location from firing rates alone, which can suffer when grid cells exhibit stable differences in in-field firing rate (Boccara, Nardin, Stella, O'Neill, & Csicsvari, 2019;Butler, Hardcastle, & Giocomo, 2019;Ismakov, Barak, Jeffery, & Derdikman, 2017). Finally, we describe analytic techniques that can be used to identify phase coding in the absence of high amplitude and continuous frequency oscillation. In particular, we show that dividing the receptive field or spike train into a small number of discrete regions and then extracting the mean phase in each region after filtering the LFP signal in a broad frequency range allows the robust identification of phase coding. Indeed, a similar approach has recently been used to identify grid and place cell phase precession against highly variable low frequency oscillations in the flying bat (Eliav et al., 2018). In summary, we demonstrate that phase coding in grid cells is not contingent on sustained rhythmicity in neural activity and offers several computational advantages over a pure rate code.

| Grid cell rate code
We simulate a total of N = 200 grid cells divided equally among m = 5 modules with a minimum scale of s 5 = 30 cm and a ratio of 1.4 between subsequent scales s m . We make use of a phenomenological model of grid cell firing in which the relative influence of rate and phase coding can be independently modulated (following Chadwick, van Rossum, & Nolan, 2015). Specifically, the firing rate of each cell r g (t) at time t is the product of a rate code r x ð x ! ) that is dictated by the agent's location x ! t ð Þ = x, y f g; and a phase code r ϕ (ϕ, θ) that is dictated by the location-dependent preferred firing phase ϕ x ! and LFP phase θ(t): The grid cell rate code r x ð x ! ) is a Gaussian function of the distance d between the agent and the center of the closest grid node x c ! = x c ,y c f g, with field width governed by a constant σ m = s m /10 and the maximum in-field firing rate given by a constant r c : In simulations of movement on a linear track, the peak locations of each grid firing field x c ! are uniformly distributed and repeat with a period equivalent to the grid scale. In simulations of movement in a two-dimensional (2D) environment, the peak locations of each grid firing field x c ! are uniformly distributed and repeat at the vertices of a rhombus with the length of each side equal to the grid scale and an acute angle of 60 .
In simulations of a uniform grid firing pattern, the value of maximum in-field firing rates r c = 1; while in simulations where firing rates vary between grid fields (as observed experimentally; Ismakov et al., 2017), the specific value for each field is drawn at random from a normal distribution with a mean and variance of one (rectified to prevent negative firing rates).

| Grid cell phase code
The grid cell phase code r ϕ (ϕ, θ) is modelled as a circular Gaussian function of the difference between instantaneous LFP phase θ(t) and the location-dependent preferred firing phase of each cell ϕ x ! , with the influence of phase coding over the spike train governed by a constant k= 1.5: In simulations of grid cells that exhibit phase precession, the location dependent preferred firing phase of each cell ϕ x ! is proportional to the linear distance from the current location to the centre of the closest grid node projected onto the direction of travel wherev is a unit vector in the direction of velocity v ! (following Burgess et al., 1994;Jeewajee, Barry, O'Keefe, & Burgess, 2008a): Conversely, in simulations of grid cells that exhibit phase locking, the preferred firing phase ϕ = π at all locations within the environment.
The overall activity of each grid cell R(t) is determined by the product of the firing rate r g (t); the instantaneous baseline frequency f (t), to normalize the number of spikes in each oscillatory cycle; and instantaneous running speed v ! multiplied by a constant m v = 0.16 cm −1 , to account for the experimentally observed increase in firing rate with running speed (Sargolini et al., 2006): Finally, spike trains for each grid cell are generated by an inhomogeneous Poisson process with rate r tot (t), which is determined by R(t) and normalized to ensure that each cell has a mean firing rate of r = 2 Hz across the duration of the simulation T:

| LFP signal
In initial simulations that aim to replicate the properties of rodent entorhinal grid cells, the LFP signal is a sinusoid with frequency f= 8 Hz. In further simulations that aim to examine the properties of putative human grid cells, the LFP signal is taken from depth electrodes in the hippocampi of pre-surgical epilepsy patients performing a spatial memory task at the National Hospital for Neurology and Neurosurgery, London, recorded at a sample rate of 512 Hz (see Bush et al., 2017 for further details). In this case, we first examine post-implantation computed tomography scans co-registered with pre-implantation magnetic resonance images to identify candidate electrode contacts. We choose those which are unequivocally located in the body of the hippocampus, as the fidelity of post-implantation images is not sufficient to confidently resolve medial entorhinal cortex. However, rodent data suggest that movement related oscillatory activity is highly coherent between these regions (Mizuseki, Sirota, Pastalkova, & Buzsaki, 2009).
Next, the LFP signal from hippocampal electrode contacts is band-pass filtered in the 2-20 Hz range using a zero-phase, secondorder Butterworth filter. The specific choice of filter band is arbitrary, as any LFP signal with broadband low frequency power from which the phases of multiple spikes can be extracted after filtering is sufficient (i.e., giving an upper limit of several hundred Hz, assuming spikes have duration of~1 ms). Next, the phase of the filtered LFP signal at each time step θ(t) is extracted using the Hilbert transform. Instantaneous frequency f(t) is then computed from the phase advance between each pair of adjacent samples and smoothed with a 50 ms box-car filter. Finally, dynamic phase and frequency information are down-sampled using linear interpolation to match the simulation time step of 5 ms (i.e., 200 Hz), and time windows that match the duration of tracking data for each run T are selected at random from randomly chosen electrode contacts. Importantly, there is no correspondence between these tracking data and the human LFP data used in each simulation, although both location (derived from the former) and LFP phase (derived from the latter) jointly determine the grid cell phase code according to Equations (3) and (4).

| Grid cell analysis
We restrict all grid cell analyses to periods of movement, defined as time bins where running speed v ! ≥ 5 cm s −1 . This is intended to match standard hippocampal electrophysiology analysis protocols, which typically exclude data from low running speeds because place and grid cells exhibit non-local coding during periods of immobility (Olafsdottir, Bush, & Barry, 2018).
First, we examine the grid cell firing rate code for location. For both 1D and 2D environments, we compute the mean firing rate in 2-cm sided bins and then smooth with a five bin boxcar kernel. Grid fields are subsequently defined as at least 5/10 contiguous bins in 1D/2D environments where firing rates are greater than 10% of the peak firing rate across the entire trial. For 2D environments, we quantify the heterogeneity of peak in-field firing rates for each cell using the coefficient of variation, which is equal to the standard deviation of peak firing rates across fields divided by their mean (Ismakov et al., 2017). In addition, rate maps are used to generate spatial autocorrelations from which gridness scores and grid scale can be estimated, as described previously (Sargolini et al., 2006). To establish whether a rate map shows significant six-fold symmetry, the true gridness is compared with the 99th percentile of a surrogate distribution generated by rotating the spike train relative to the tracking data by a time shift sampled from a random uniform distribution in the range  Eliav et al., 2018;Royer, Sirota, Patel, & Buzsáki, 2010). This is achieved by setting the value of the autocorrelation at zero lag to the maximum value across all time lags t and then fitting the corrected autocorrelation with a function A(t): Fitted parameter values are restricted to the following ranges: Nonlinear least squares fitting is performed 500 times using the Matlab "fit" function with different initial values drawn randomly from a uniform distribution within the ranges specified above, and the estimated values of each parameter taken from the fit with the greatest R 2 value compared to the true autocorrelation. The oscillation index is then computed as a/max(A(t)).
To compute the instantaneous frequency difference between simulated grid cell spikes trains and the LFP, we estimate the oscillatory phase of the mean normalized grid cell phase code r ϕ (ϕ, θ) at each time step using the Hilbert transform, compute instantaneous frequency as the change in phase between time steps, and then compare this to the LFP frequency within the same time bin. Finally, phase precession is quantified as the circular-linear correlation between firing phase and distance travelled through the grid field (Jeewajee, Barry, et al., 2008a;Kempter, Leibold, Buzsáki, Diba, & Schmidt, 2012). The significance of this correlation is established by randomly shuffling the firing phase values relative to the distance values 1,000 times, re-computing the value of the circular-linear correlation each time and then comparing the true correlation value to this surrogate distribution.

| Decoding
Having simulated grid cell population activity, we then attempt to reconstruct the running speed, location and movement direction of our simulated agent, as well as an arbitrary fourth variable, from grid cell population firing patterns within each oscillatory cycle. To do so, we first compute the total number of spikes fired by each cell k i ! in oscillatory cycle i (i.e., the N dimensional population vector); as well as the average location x i ! and running speed v i ! of the simulated agent in that cycle.
We then decode running speed as follows: we estimate the slope and intercept of the relationship between the total number of spikes fired (i.e., the sum of population vector k i ! across all cells) and running speed in each cycle, using data from alternate cycles to avoid overfitting. We predict running speed in the remaining oscillatory cycles based on the total number of spikes fired by the grid cell population, and quantify decoding accuracy by comparing predicted and actual running speed values.
We then decode location, irrespective of firing phase, for all cycles with an average running speed ≥5 cm s −1 using maximum likelihood estimation (following Mathis et al., 2012). To do so, we compare the population vector k i ! from cycle i with the expected firing rate in each location r x x ! , independent of phase coding, produced by averaging the grid cell rate for location r x x ! across 2 cm sided bins covering the entire environment, multiplied by the average cycle duration T cyc . The maximum likelihood estimate of the agent's location given the population vector from that cycle x ! MLE k i ! is then given by: We decode location, incorporating phase information, as follows: we divide each oscillatory cycle i into n φ = 5 phase bins with approximately equal spike counts (again, excluding cycles with an average running speed <5 cm s −1 ) and compute the N × n φ dimensional population vector k i,φ ! (i.e., the total number of spikes fired by each cell in each phase bin of each oscillatory cycle). We then compute the expected number of spikes in each phase bin of each cycle based on the combined rate and phase code for location r g multiplied by the average phase bin duration T φ . As a control, to quantify the specific contribution of phase information, we also compute the expected number of spikes in each phase bin of each cycle based on the rate code for location alone r x multiplied by the average phase bin duration T φ . The maximum likelihood estimate of the oscillatory cycle that produced the observed population vector i g k i,φ ! based on a combined rate and phase code is then given by: In addition, the maximum likelihood estimate of the oscillatory cycle that produced the observed population vector i x k i,φ ! based on a pure rate code is given by: In both cases, the decoded location is taken as the average location in the decoded oscillatory cycle. Importantly, in the presence of stable variation in peak in-field firing rates, we can construct a "naïve" decoder by setting r c = 1 when computing r g and r x above; and an "informed" decoder by setting r c to the values used to generate grid firing patterns for each cell.
Finally, we decode movement direction as follows: we use maximum likelihood estimation to decode location independently within each phase bin of each oscillatory cycle as described above (i.e., using Equation (8), but replacing the N dimensional population vector k i ! from cycle i with the N dimensional population vector k i,p ! from phase bin p in cycle i). This generates a sequence of decoded locations within each oscillatory cycle. We compute the slope of x MLE, φ and y MLE, φ across phase bins using linear regression, and estimate movement direction as the inverse tangent of those slopes. The decoded movement direction can be compared with the true movement direction computed from the actual location in each using the same method.

| Empirical tests of phase coding
We also use our simulated data to examine how phase coding might be identified in typical empirical data sets consisting of spike times and an aperiodic LFP signal. To do so, we examine spike triggered LFP averages at different locations within the firing field, consistent with previous studies (Eliav et al., 2018). Importantly, however, we wish to eliminate any differences in amplitude between cycles, which can skew spike triggered LFP averages and thus reduce the sensitivity of this approach. We therefore generate a synthetic LFP signal s(t) with constant amplitude from the phase of the filtered LFP signal at each time step θ(t): We then compute the average of this synthetic LFP signal in three equally sized regions of each firing field for each cell, extract the circular mean phase at the time of firing across fields and cells for analysis purposes, and the average synthetic signal across fields and cells for illustration purposes.

| Approximating the LFP signal
Finally, we examine whether population spiking activity could provide a good substitute for the LFP baseline signals used in these simulations. To do so, we sum the firing rates of all simulated grid cells in each temporal bin and then band-pass filter that multi-unit activity using the same 2-20 Hz second-order Butterworth filter applied to the true LFP signal (see above). In this case, however, we make use of a causal filter to approximate a biologically realistic leaky integration process (i.e., which can only depend on past, and not future, inputs).
We then correlate filtered multi-unit activity with the synthetic LFP signal s(t), to avoid confounds arising from dynamic changes in amplitude.

| Data availability
Data available on request from the authors. Code for all simulations, along with sample LFP data, is available at http://modeldb.yale.edu/ 261878.

| Grid cell phase coding in the absence of rhythmicity
In contrast to the rodent, prominent and sustained oscillatory activity in LFP recordings from the human brain are rare. This has led to the suggestion that phase coding cannot play a role in human cognition, despite the numerous theoretical advantages afforded by such a scheme. Here, using intracranial LFP recordings from the human hippocampus, we aim to demonstrate that a robust phase code can be maintained in the absence of a prominent baseline oscillation with relatively fixed frequency; and to examine the functional advantages offered by phase coding over and above a pure rate code. As a model system, we make use of grid cells, which are characterized by a regular triangular array of spatial firing fields (Hafting, Fyhn, Molden, Moser, & Moser, 2005). In the rodent, grid cell firing is also modulated by the movement related 6-12 Hz LFP theta oscillation and a significant proportion of grid cells exhibit theta phase precession-firing at successively earlier phases of each theta cycle as their firing field is traversed (Climer et al., 2013;Hafting et al., 2008;Jeewajee et al., 2014).
In humans and flying bats, however, grid cells appear to exist in the absence of a prominent LFP oscillation (Eliav et al., 2018;Jacobs et al., 2013;Yartsev et al., 2011).
We begin by simulating the firing patterns of rodent grid cells using a phenomenological model in which the influence of location on both firing rate and preferred firing phase can easily be modulated (Chadwick et al., 2015; see Section 2). An examination of firing rates confirms that simulated grid cells show periodic spatial firing patterns in both 1D ( Figure 1a) and 2D (Figure 1b-d) environments. This periodicity can be quantified using the gridness metric, which is significant for all cells in the 2D environment, but lower for larger scale grids that The grid cell rate code. As observed in vivo, simulated grid cells exhibit a firing rate code for location that corresponds to periodic spatial firing fields both (a) on a linear track; and (b) in the open field. This is confirmed by an inspection of the (c) firing rate map (peak rate inset); and (d) spatial autocorrelation (gridness score inset). (e) Distribution of gridness scores across the population, which are significant for all cells (median value inset). Lower gridness scores are produced by larger scale grids, which typically exhibit only two firing fields within the environment. (f) Finally, grid cells show substantial variation of peak in-field firing rates, quantified using the coefficient of variation (CV) for each cell (median value inset). Panel a shows data from a single simulated grid cell, averaged over 20 independent simulations; panels b-d show data for a single simulated grid cell from a single simulation; panels e and f for the entire grid cell population from a single simulation [Color figure can be viewed at wileyonlinelibrary.com] exhibit fewer fields within the arena (Figure 1e). Interestingly, simulated grid cells also exhibit substantial variation of peak in-field firing rates, despite these simulations being intended to produce a uniform grid firing pattern (Figure 1f). This is likely to result from Poisson firing statistics and non-uniform coverage of the environment, consistent with previous theoretical studies (Ismakov et al., 2017).
As in the rodent entorhinal cortex, grid cell firing in these initial  (Eliav et al., 2018). Importantly, however, the burst firing frequency of grid cells is slightly higher than LFP theta, and this difference increases with running speed at a rate that is inversely proportional to the measured grid scale (Geisler et al., 2007(Geisler et al., , 2010; Figure 2f). As a result, the average frequency difference between the burst firing of grid cells and the LFP across the entire simulation is inversely proportional to grid scale (Figure 2g; Jeewajee, Barry, et al., 2008a). Hence, grid cells across the population exhibit theta phase precession-firing progressively earlier in each theta cycle as their firing fields are traversed (Figure 2h,i).
Next, we simulate the firing pattern of grid cells recorded in the human or flying bat brain where prominent, sustained oscillations are rare (Eliav et al., 2018;Jacobs et al., 2013;Yartsev et al., 2011). To do so, we make use of LFP recordings extracted from depth electrodes in the hippocampal formation of pre-surgical epilepsy patients performing a spatial memory task in place of a constant 8 Hz theta oscillation (Bush et al., 2017). In this case, LFP frequency varies dynamically over a broad range, such that no clear peak is observed in either the power spectrum (Figure 3a)  In each of the simulations described above, all grid cells exhibit phase precession-that is, their preferred firing phase is dictated by progress through the firing field (see Section 2). However, in empirical data from both rodents (Climer et al., 2013;Hafting et al., 2008;Jeewajee et al., 2014;Sargolini et al., 2006) and flying bats (Eliav et al., 2018), a significant proportion of grid cells also exhibit phase lockingconsistently firing at the same phase of low frequency fluctuations in the LFP. Hence, we next simulate grid cell firing that is phase locked to the highly variable human intracranial LFP signal, as these data will later serve as a useful control to establish the specific functional contribution offered by phase coding. In these simulations, gridness scores remain In sum, these results demonstrate that both the phase precession and phase locking of grid cell firing could occur in the absence of any apparent rhythmicity in either the LFP or spike train, given that phase is independent of frequency (Blair et al., 2013;Orchard, 2015). This is consistent with recent empirical data from bats, where place and grid cells show either phase locking or phase precession relative to an LFP signal with a frequency that varies dynamically over a wide range (Eliav et al., 2018). Importantly, this also demonstrates that LFP dynamics in the human hippocampus do not preclude the existence of robust phase coding (see Section 4 for potential mechanisms). Nonetheless, it should be emphasized that none of the results presented here are contingent on specific features of those LFP data. Any baseline signal with broadband low frequency power from which phase information can be extracted after filtering (here, in the 2-20 Hz range) would produce the same qualitative output, regardless of whether that signal exhibited rhythmicity within a specific, narrow frequency range. Indeed, multi-unit activity in these simulations provides a good substitute for the LFP signal (see Section 2 and Figure 7d,e). It is also important to note that there is no correspondence between the rodent behavioral trajectories and human LFP data used in these simulations, although both location (derived from the former) and LFP phase (derived from the latter) jointly determine the grid cell phase code for location. Next, we examine the potential functional contribution made by that phase code in grid cells; and how this might be affected by the absence of a prominent baseline oscillation with approximately constant frequency.

| Grid cells multiplex spatial information in firing rate and phase
To address the potential contribution made by phase coding to cognition, we next examine what information can be decoded from grid cell population activity within each oscillatory cycle and how this is affected by the presence or absence of a phase code for location within the firing field. Several previous theoretical studies have established that, while the firing pattern of a single grid cell or module of grids cells is inherently ambiguous about an animal's location, grid F I G U R E 3 The grid cell temporal code in the absence of rhythmicity. In these simulations, we make use of a local field potential (LFP) signal recorded from the human hippocampal formation with a frequency that varies dynamically over a broad range, such that no clear peak is observed in either the: (a) LFP power spectrum (6-10 Hz theta band marked in gray); or (b) spike train temporal autocorrelation (oscillation index inset). (c) Distribution of oscillation indices, which are significantly lower than simulations with constant theta rhythmicity (median value inset).
(d) Distribution of gridness scores across the population, which are significant for 99.5% of cells and do not differ from simulations with constant theta rhythmicity (median value inset). (e) Distribution of peak in-field firing rate coefficient of variation (CV) values, which do not differ from simulations with constant theta rhythmicity (median value inset). (f) Distribution of firing phase resultant vector lengths across the population, with 99% of cells showing significant phase modulation (median value inset). (g) Correlation between distance travelled through the grid field and firing phase for a typical grid cell, illustrating phase precession (circular-linear correlation coefficient inset, and line of best fit plotted in red).
(h) Distribution of circular-linear correlation coefficients across the population, which are significant for all cells and do not differ significantly from those with constant theta rhythmicity (median value inset). (i) The intrinsic firing frequency of individual grid cells continues to exceed the LFP frequency by an amount that is inversely proportional to their grid scale, despite that frequency varying dynamically over a wide range. Panel a shows data from a single simulation; panels b and g for a single simulated grid cell; panels c-f, h, and i for the entire grid cell population [Color figure can be viewed at wileyonlinelibrary.com] cell population activity across modules with different scales can offer a robust and accurate code for location and navigation over a large range (Bush, Barry, Manson, & Burgess, 2015;Fiete et al., 2008;Mathis et al., 2012;Stemmler, Mathis, & Herz, 2015). Consistent with those studies, we find that even the relatively small population examined here (five modules, containing a total of 200 grid cells) can encode location with a median error of <2 cm in both 1D ( Figure 5a) and 2D (Figure 5b) environments. In the larger (~50 m) 1D environment, however, the grid cell rate code is prone to occasional catastrophic errors (defined here as ≥50 cm), consistent with previous studies (Fiete et al., 2008;Mathis et al., 2012;Towse, Barry, Bush, & Burgess, 2014). These large errors reflect the upper capacity limit of the grid cell network, which is generally dictated by the lowest common multiple of grid scales (Bush et al., 2015;Fiete et al., 2008; but see Mathis et al., 2012;Stemmler et al., 2015).
In addition to encoding location, empirical data indicate that the mean firing rate of grid cells-like those of other principal cells in the hippocampal formation-varies with running speed (Sargolini et al., 2006). The cells simulated here also exhibit this property (see Section 2), such that we can estimate the agent's running speed from the total number of spikes fired by the grid cell population within each oscillatory cycle (as demonstrated previously for speed cells; Gois & Tort, 2018;Kropff, Carmichael, Moser, & Moser, 2015). This allows us to decode running speed with an accuracy of ≤5 cm s −1 in 95% of cycles ( Figure 5c). Hence, these results demonstrate that both location and movement speed can be accurately decoded from the firing rate of a small number of grid cells in each oscillatory cycle and, importantly, that the modulation of firing rates according to running speed does not compromise location decoding.
Next, we turn to the nature of information encoded by the grid cell phase code, which has previously received relatively little attention. In place cells, it has been demonstrated that theta phase precession causes sweeps of activity corresponding to the animal's current trajectory to occur within each oscillatory cycle (Burgess et al., 1994;Johnson & Redish, 2007;Skaggs et al., 1996; Figure 5d). It has been suggested that movement direction might therefore be encoded by the relative phase of place or grid cell firing across the population (Zutshi et al., 2017), although this has yet to be demonstrated. Consistent with this hypothesis, we demonstrate that-by fitting a linear trajectory to the sequence of locations decoded independently from five phase bins within each oscillatory cycle-we can estimate movement direction with an accuracy of ≤30 in 76% of cycles ( Figure 5e). In contrast, the decoding of movement direction is significantly less accurate in simulations without phase coding, where an accuracy of ≤30 is achieved in only 30% of cycles (Figure 5f), although still possible in some cases due to the small but significant movement made within each oscillatory period. Specifically, the distribution of movement direction decoding errors is non-uniform with a mean of zero both with (V = 4,745, p < 0.001) and without (V = 1,742, p < 0.001) phase coding, although the variance is significantly greater in the latter case (k = 18.8 × 10 6 , p < 0.001). This demonstrates that the grid cell phase code for location contributes to the accurate encoding of movement direction in population activity within each oscillatory cycle.
F I G U R E 5 Decoding movement trajectories from grid cell population activity in each oscillatory cycle. Location in both (a) one-dimensional (1D) and (b) two-dimensional (2D) environments can be decoded from population firing rates, using the grid cell rate function (frequency of catastrophic errors inset). (c) In addition, running speed can be decoded from the total number of spikes fired by the grid cell population (standard deviation inset). (d) Phase precession across the grid cell population generates "theta sequences" of activity within each oscillatory cycle that correspond to the current movement trajectory (color axis indicates spike density). (e) This allows movement direction to be estimated from the sequence of locations decoded across phase bins (resultant vector length inset). (f) Conversely, in the absence of a phase code for location, movement direction decoding is much less accurate (resultant vector length inset). (g) Location can also be decoded from population firing rates, using the history of firing rates across oscillatory cycles in that simulation (frequency of catastrophic errors inset). (h) However, decoding location using both population firing rate and phase information is more accurate, as illustrated by a significant reduction in the rate of catastrophic errors (frequency of catastrophic errors inset). Panels a, b, g, and h are averaged over data from 20 independent simulations [Color figure can be viewed at wileyonlinelibrary.com] Importantly, the decoding of movement direction from grid cell population activity does not require an explicit measure of grid cell firing phase-downstream neurons could simply read out a location estimate from simultaneously active grid cells in different phase bins by coincidence detection. Movement direction is subsequently estimated from the sequence of decoded locations across the oscillatory cycle, and each sequence of decoded locations is punctuated by a period of relatively low population firing rates between oscillatory cycles. Similarly, accurate decoding of movement direction does not require the phase of firing in individual grid cells to be particularly precise-only that firing phase can be measured. In these simulations, the combination of Poisson firing statistics and a relatively broad preferred firing phase distribution (e.g., with a circular standard deviation of~0.9 rad, see Equation (3)) produce a significant level of phase noise in simulated grid cells. Moreover, the results described above are generated by independently decoding location from grid cell population activity in five phase bins, corresponding to a phase resolution of~72 or 2π/5.
Elsewhere, previous studies have demonstrated that incorporating phase information can improve the accuracy with which location is decoded from place cell firing rates (Jensen & Lisman, 2000). To establish if this is also true for grid cell firing patterns, we compare the accuracy with which location can be decoded in the larger (~50 m) 1D environment. We find that incorporating the phase of firing significantly improves location decoding accuracy, compared to decoding using firing rates alone (Figure 5g,h). Specifically, incorporating phase information significantly reduces the incidence of catastrophic location decoding errors (t[38] = −12.0, p < 0.001, d = 3.79). Finally, we note that the frequency of the baseline oscillation in these simulations varies independently of grid cell firing rate and phase, and could therefore be used to encode an additional variable. In the rodent, for example, it has been demonstrated that LFP theta frequency encodes running speed information, like neural firing rates, and possibly driven by the same medial septal glutamatergic inputs (Fuhrmann et al., 2015;Hinman, Brandon, Climer, Chapman, & Hasselmo, 2016;Wells et al., 2013). Similarly, in the flying bat, it has been demonstrated that LFP frequency is higher during faster movements (Eliav et al., 2018).
Although baseline frequency is not modulated by running speed in these simulations, we find that it is possible to accurately decode an arbitrary fourth variable from the duration of each oscillatory cycle by linear regression, with an error of ≤5% in 91% of cycles (data not shown).
In sum, these results demonstrate that the rate and phase code for location provided by grid cells can efficiently and accurately multiplex information about the complete movement trajectory of an agent, including location, running speed, and movement direction, in addition to a fourth variable provided by the oscillatory period. In the absence of phase coding, movement direction can no longer be accurately decoded from grid cell activity, and the accuracy of location decoding is reduced. Hence, phase coding in grid cells makes a specific contribution to the encoding of movement information above and beyond that provided by an examination of firing rates alone. Crucially, each of these results could equally arise from decoding place cell firing patterns, which exhibit an equivalent rate and phase code for location. Our decision to simulate grid cell firing patterns was driven solely by computational efficiency, as fewer grid cells are required to accurately encode location (Fiete et al., 2008) and inherently provide a greater number of runs through the firing field in each simulation.

| Phase coding with variable in-field firing rates
Interestingly, it has been demonstrated that grid cells exhibit stable differences in peak in-field firing rates across familiar environments (Ismakov et al., 2017), with grid fields that are closer to persistently rewarded locations generally being more active (Boccara et al., 2019;Butler et al., 2019). It has been suggested that this allows salient locations to be encoded more accurately by the grid cell population, although empirical proof is lacking. We therefore chose to examine the impact of firing rate variability on the accuracy with which location can be decoded from the grid cell population by simulating grid cells with stable differences in peak in-field firing rate that approximate the variability observed in vivo (Ismakov et al., 2017). As In this case, firing rate variability introduces ambiguity into the grid cell rate code for space, as low firing rates may indicate either the periphery of a field with high peak firing rate or the center of a field with low peak firing rate.
Intuitively, this ambiguity might be resolved by incorporating phase information, which specifically indicates the location of an agent within the firing field. Indeed, we find that decoding location using a uniform firing rate function alongside information about the phase of firing produces fewer catastrophic errors than the naïve rate decoder (t[18] = −19.7, p < 0.001, d = 6.22; Figure 6i, compare with F I G U R E 6 Phase contributes to accurate location decoding with variable in-field firing rates. (a) Sample rate map from a simulated cell with stable differences in peak in-field firing rate (overall peak firing rate inset). (b) Corresponding spatial autocorrelation (gridness score inset).
(c) Distribution of gridness scores across the population, which are significant for 94% of cells but significantly lower than simulations with a uniform grid firing pattern (median gridness score inset); (d) Distribution of peak in-field firing rate coefficient of variation (CV) values, which are significantly higher than simulations with a uniform grid firing pattern (median value inset); (e) Distribution of firing phase resultant vector lengths, which are significant for all cells (median value inset); and (f) Distribution of location versus phase correlation coefficients across the population, which are significant for all cells (median value inset). (g) Location decoding with an "informed" decoder. Accuracy is significantly improved when in-field firing rates are variable and the decoder is informed about the expected firing rate in each field (compare with Figure 5g, frequency of catastrophic errors inset). (h) Location decoding with a "naïve" decoder. Accuracy is significantly worse when the decoder is naïve to the expected firing rate in each field (compare with Figure 5g, frequency of catastrophic errors inset). (i) Location decoding with a "naïve" decoder and phase information. Providing the naïve decoder with information about firing phase reduces ambiguity without the need for environment specific learning (frequency of catastrophic errors inset). Panels a and b show data for a single simulated grid cell; panels c-f for the entire grid cell population; and panels g, h, and i are averaged over data from 20 independent simulations [Color figure can be viewed at wileyonlinelibrary.com] Figure 6h), consistent with the results described above (Figure 5g,h).
Although the incidence of catastrophic errors remains slightly but significantly greater than that produced by an informed decoder (t [38] = 11.7, p < 0.001, d = 3.69; compare with Figure 6g), the agent does not need to learn the configuration of peak firing rates in each novel environment. This is consistent with the notion that grid cells support the generalization of structural regularities between contexts (Behrens et al., 2018). In sum, these results indicate that variable infield firing rates only contribute to more accurate encoding of location by grid cells if the decoding apparatus is informed of the expected variability; and that providing the decoding apparatus with access to phase information can ameliorate this issue by resolving location within the firing field. This highlights another potential contribution of phase coding to cognition: disambiguating stimuli (i.e., locations) that produce similar firing rates.

| Empirical tests for phase coding in the absence of rhythmicity
Finally, we consider the problem of how phase coding can be identified in empirical data when sustained rhythmicity is absent. One F I G U R E 7 Empirical tests for phase coding in the absence of rhythmicity. (a) Spike triggered average of the local field potential (LFP) signal in early, middle and late parts of the firing field (one-dimensional [1D] environment). (b) Spike triggered average of the same LFP signal, reconstructed from phase information after a broadband filter is applied, in early, middle and late parts of the firing field (1D environment). Peak deflection of the LFP signal progresses from a negative to positive time lag with progress through the receptive field, consistent with the presence of a phase code for location within the field. (c) Circular mean phase at the time of firing in early, middle and late parts of the firing field, averaged across all cells. Note that this is equivalent to the phase precession plots shown in Figures 2h and 3g (Eliav et al., 2018). Complications arise, however, from the fact that LFP amplitude may vary independently of frequency, such that the LFP signal associated with spikes fired during periods of relatively high LFP amplitude can dominate spike-triggered averages (Figure 7a). To avoid this potential confound, the baseline signal can be reconstructed from the cosine of its phase, which is extracted from the analytic signal after a broadband filter has been applied-effectively orthogonalizing phase information from variations in amplitude.
By using this approach on our simulated data, we find that spiketriggered LFP averages clearly differ between early, middle and late parts of the grid firing field when simulated cells exhibit phase precession, but not when they exhibit phase locking ( Figure 7b). Indeed, the circular mean phase at the time of spiking changes significantly across early, middle and late parts of the receptive field (Figure 7c), consistent with the phase precession analyses described above (e.g., Figures 2h and 3g). We emphasize that this approach is not specific to spatially modulated firing patterns-similar plots could be generated for different stimulus intensities, or positions within a spike train, to assay the presence of phase coding in a diverse range of cortical regions (e.g., Aghajan et al., 2015).

| DISCUSSION
The simulations presented here support two important conclusions regarding the potential role of phase coding in human cognition. First, that phase coding does not necessarily rely on a prominent baseline oscillation with relatively invariant frequency, and thus that the absence of such a signal in typical LFP recordings from the human brain does not preclude a role for phase coding. This is consistent with recent empirical data from the flying bat, where a phase code for location is observed in place and grid cells of the hippocampal formation, despite the absence of any sustained oscillatory activity in a fixed frequency band (Eliav et al., 2018). Second, that phase coding allows for both the multiplexing of additional information in spike trains and the disambiguation of stimuli that generate similar firing rates, and thus offers computational advantages above and beyond those provided by firing rates alone. This is consistent with previous theoretical studies (Fries et al., 2007;Jensen et al., 2014;Panzeri et al., 2010;Thorpe et al., 2001) and empirical data (Kayser et al., 2009;Montemurro et al., 2008;O'Keefe & Burgess, 2005;Siegel et al., 2009;Turesson et al., 2012;Zuo et al., 2015), which demonstrate that phase modulation permits greater information content.
In the context of spatial cognition, we have demonstrated that phase coding within each oscillatory cycle encodes movement direction, consistent with previous suggestions (Zutshi et al., 2017).
Movement direction cannot be reliably recovered from firing rates alone, although these can robustly encode both location and movement speed (Fiete et al., 2008;Mathis et al., 2012 (Raudies, Brandon, Chapman, & Hasselmo, 2015). Hence, the firing phase of grid and place cells in the hippocampal formation is the only robust correlate of movement direction identified thus far, raising the question of which upstream circuits provide the requisite input.
Although these simulations demonstrate that prominent and sustained rhythmicity is not necessary to support phase coding, they do not address the question of why an arrhythmic state might be preferred in the human brain (Bush & Burgess, 2019). We have shown that it is possible to encode an additional variable in the ongoing frequency of the baseline signal, although this does not necessitate variability over a wide frequency range. Indeed, running speed is encoded in the frequency of rodent hippocampal theta oscillations, which typically vary by <1 Hz (Jeewajee, Barry, et al., 2008a). An alternative possibility is that changes in baseline frequency are used to switch between encoding and retrieval modes of operation, modulating synaptic plasticity by adjusting the temporal interval between firing in connected cells without affecting their relative phase. Indeed, it has been demonstrated that environmental novelty-which might be associated with enhanced learning-reduces theta frequency in the hippocampus (Jeewajee, Lever, Burton, O'Keefe, & Burgess, 2008b) by reducing the slope of the theta frequency versus running speed relationship (Wells et al., 2013). In addition, it has been proposed that encoding and retrieval processes may be segregated by oscillatory phase (Hasselmo, Bodelon, & Wyble, 2002), and this is supported by recent empirical work showing that the preferred theta firing phase of hippocampal place cells also changes during novelty (Douchamps, Jeewajee, Blundell, Burgess, & Lever, 2013).
Phase coding may also offer functional advantages for cognition beyond those considered here. In the context of the rodent hippocampal formation, for example, it has recently been demonstrated that the offline "replay" of behavioral trajectories in place cell activity relies on robust theta sequences during learning (Drieu, Todorova, & Zugaro, 2018). Interestingly, the emergence of robust theta sequences during development also coincides with the disappearance of infantile amnesia (Farooq & Dragoi, 2019;Muessig, Lasek, Varsavsky, Cacucci, & Wills, 2019;Travaglia, Bisaz, Sweet, Blitzer, & Alberini, 2016). These results suggest that the temporal organization of hippocampal spiking activity during active navigation contributes to robust memory encoding. In addition, theta sequences might be useful for the prospective evaluation of upcoming locations or trajectories during active navigation (Bicanski & Burgess, 2018;Bush et al., 2015;Johnson & Redish, 2007). Hence, the phase coding of sensory information during active experience is also likely to contribute to the flexible planning of subsequent behavior.
Finally, although these results demonstrate that it is theoretically possible to maintain a robust phase code in the absence of a baseline oscillation of approximately constant frequency, it is also important to consider how such phase coding might be practically established in real neural circuits. In these simulations, we have assumed that the intrinsic firing frequency of grid cells exceeds that of the baseline oscillation by a value that is proportional to running speed and inversely proportional to grid scale (Figures 2g and 3i). Because phase is the time integral of frequency, and distance the time integral of running speed, this ensures that firing phase encodes displacement through the firing field. Hence, any mechanism that can dynamically maintain the intrinsic firing frequency of a cell above baseline by an amount that is proportional to running speed through the firing field is sufficient to produce robust phase coding. The most parsimonious explanation for this precise frequency difference between baseline and active oscillatory signals being maintained over long periods is that the baseline frequency (assumed here to be represented by the LFP) reflects the average frequency of velocity controlled oscillatory inputs with different preferred firing directions (Burgess, 2008;Burgess, Recce, & O'Keefe, 1993;Geisler et al., 2007Geisler et al., , 2010Hasselmo, 2008;Welday, Shlifer, Bloom, Zhang, & Blair, 2011). Under such a scheme, the LFP signal effectively reflects input from a population of cells whose intrinsic firing frequency varies as a cosine function of movement direction.
Importantly, this proposed mechanism for phase coding need not be restricted to spatial cognition (e.g., Terada, Sakurai, Nakahara, & Fujisawa, 2017)-any variable can be encoded in the phase of neural firing if the difference between intrinsic firing frequency and LFP baseline frequency is proportional to the time derivative of that variable. This raises the possibility that differential frequencies in active neurons and the LFP represent a domain general mechanism in support of phase coding. Crucially, as described above, this frequency difference can be established and dynamically maintained if the LFP signal simply corresponds to the smoothed average population activity (Burgess et al., 1993;Geisler et al., 2007Geisler et al., , 2010. This is consistent with the hypothesis that extracellular voltage recordings from a specific region (i.e., the LFP) primarily reflect net synaptic input to that region, if we assume that population activity is a linear function of that input (Buzsaki, Anastassiou, & Koch, 2012). Indeed, in these simulations, filtering the multi-unit activity of all grid cells in the same 2-20 Hz low frequency band as the LFP signal consistently provides a good substitute for that signal (Figure 7d,e, see Section 2). In sum, we have demonstrated that phase coding in grid cells is not contingent on sustained rhythmicity in neural activity and offers several computational advantages over a pure rate code. Hence, it seems likely that the arrhythmic activity patterns observed in the human brain do not preclude a role for phase coding, and that such a code could make a unique contribution to higher cognitive functions.