Boundary tuning while changing environmental geometry

To test the hypothesis that the recorded cells were tuned to a preferred boundary rather than to a specific place field in the environment, we conducted an additional experiment. In this control experiment, we measured the neural activity of the same cells before and after we changed the tank geometry. This was done by adding an additional horizontal half-wall (control in 8 of the 35 boundary vector cells). Examples of these cells are presented in Fig 5. Each cell was recorded for a full recording session (60 to 75 min; Fig 5A), after which the fish was blindfolded, and a Perspex shelf was inserted into the water tank (Fig 5B). Another full recording session was then initiated.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Changing environmental geometry.

(A) The first session was recorded in the main experimental water tank. (B) Before the second recording session, a shelf was inserted into the water tank to modulate the geometry of the environment. (C) Example rate map of a boundary vector cell tuned to distance from the bottom of the water tank, color coded from dark blue (zero firing rate) to dark red (maximal firing rate, indicated on the top right side of the panel). (D) Rate map over the water tank of the same cell after adding the shelf. Firing was modulated by both the bottom of the water tank and above the shelf. (E-H) Additional examples of boundary vector cells before and after the geometric change in the environment. The underlying data supporting all panels in this figure can be found in a file named Fig 5_data.mat (see Data Availability).

https://doi.org/10.1371/journal.pbio.3001747.g005

In one example, in the first recording session, the cell’s neural activity gradually decreased with distance from the bottom of the water tank (Fig 5C). After adding the shelf (gray mark in Fig 5D), the cell’s rate map (Fig 5D) shows it responded to both the bottom of the water tank and the Perspex shelf, as expected from a boundary vector cell. Additional examples showing firing patterns of boundary vector cells before and after the geometric change in the environment are presented in Figs 5E–5H and S3A–S3G.

Population analysis of boundary vector cells

The firing patterns observed in the recorded population of boundary vector cells mostly exhibited a gradual monotonous decreasing pattern with distance from a boundary found at the preferred direction of each cell (Θ_max; see Materials and methods), that is regardless of the distance of the fish from a boundary found in the orthogonal direction (Θ_max±90^o). Therefore, we calculated the receptive field size (Fig 6A) of these cells as the width of the sigmoid curve, which best fit the data (see Materials and methods and S4 Fig). The results showed that these widths distributed around 67% ± 53% (mean ± standard deviation) of the width of the water tank (Fig 6A, Units: 1 [Aq] = 70 [cm]). This fit differs from a simple border cell, which is expected to have sharp sigmoidal tuning near a specific boundary in the environment [25]. It also differs from the edge cells found previously in goldfish [14], which fired in short distance from all four walls of a shallow environment. Specifically, none of cells presented here could passed these criteria for edge encoding as in [14]. Altogether, this suggests that the cells presented here can be best described as boundary vector cells, each of which is tuned to a specific boundary direction and gradually encodes the distance from it. The Θ_max directions across the recorded population were quasi-uniformly distributed in space (Fig 6B, color-coded by the cells’ max-p values; see Materials and methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Population analysis of boundary vector cells.

(A) Distribution of receptive field widths (i.e., the width of the sigmoid curves fitted to the boundary vector cells; see Materials and methods) relative to the width of the water tank. Units: 1 [Aq] = 70 [cm]. (B) The population max-p values and preferred boundary direction Θ_max (see Materials and methods) of all 119 single units recorded from the brains of 15 fish. The dots are distributed on a radial plot where folded to the range [0, 5] and Angle = Θ_max. The color bar on the right-hand side of panel B spans the population max-p values on a logarithmic scale and shows the distribution of cells in panels B-G. (C-G) Population statistics. Y-axis in panels C-F is the rank-order max-p value of each cell (see Materials and methods), so that the cells with the strongest spatial modulation are in the top of each panel. Statistics summary of these panels is indicated in the bottom right side of each panel. For each cell, we show the following properties: (C) in-session stability correlation coefficient, (D) spatial coherence, (E) spatial information, (F) mean firing rate, and (G) mean action potential amplitude and full width at half maximum for the electrode recorded the strongest spikes. Two-sample t test results are indicated. (H) Mean firing rate in the 30% strongest bins in the rate maps of the 35 boundary vector cells in the first vs. second half of the experiments. Paired t test results are indicated. (I) Preferred boundary direction (Θ_max) of the 35 boundary vector cells vs. peak firing rate. The angle in the polar plot is Θ_max, and the radius is the peak firing rate. Circular-linear correlation coefficient (r) and p-value are indicated. The underlying data supporting all panels in this figure can be found in a file named Fig 6_data.mat (see Data Availability).

https://doi.org/10.1371/journal.pbio.3001747.g006

The characteristics of the entire recorded population are presented in Fig 6C–6G using the same color-code as in panel B. Y-axis in panels C-F is the rank-order max-p of each cell (see Materials and methods), so that the cells with the strongest spatial tuning are at the top of each panel. Rank-order correlation coefficient (r) and p-values are indicated to the bottom right of each one of these panels. As expected from space encoding cells, stronger spatial tuning (i.e., smaller max-p value) was characterized by a stronger in-session spatial-stability correlation-coefficient (Fig 6C, r = 0.67, p < 1 × 10⁻³), stronger spatial coherence (Fig 6D, r = 0.82, p < 1 × 10⁻³), and stronger spatial information (Fig 6E, r = 0.37, p < 1 × 10⁻³). In addition, the strength of the spatial modulation was also correlated with the mean firing rate during the entire recording session (Fig 6F, r = 0.45, p < 1 × 10⁻³).

Furthermore, the mean action potential properties of each cell (amplitude versus width at half maximum for the electrode that recorded the strongest spikes) are presented in Fig 6G. We used a two-sample t test and found that the ratio between the peak spike amplitude and its width can be distinguished between the 35 boundary vector cells and the rest of the recorded population (p < 1 × 10⁻³, t = 3.57 and degrees of freedom = 117).

For the 35 boundary vector cells, we also show the mean firing rate in the 30% strongest bins of each rate map (first half versus second half of each experiment; Fig 6H). A paired-sample t test (p = 0.011, t = 2.67 and degrees of freedom = 32) suggested the firing rate in the receptive field of the boundary vector cells might decrease over time. This might be related to adaptation of the fish to the environment, but further study is needed to establish this point. Last, we also show a lack of relationship between the preferred boundary direction of these cells and their peak firing rates (Fig 6I; circular-linear correlation-coefficient, r = 0.25 and p = 0.14).

The boundary vector cells’ receptive fields covered the entire water tank. A rate map overlap of the top 50 percentile bins of all boundary vector cells is presented in S3H–S3J Fig (divided into 3 different panels by the preferred direction solely for ease of visualization). The full coverage suggests that boundary vector cells might be sufficient for the fish to encode position in small environments.

We wanted to verify that the spatial firing pattern of these cells is, indeed, due to position encoding rather than other aspects of locomotion, such as swimming direction or swimming speed. For example, a gradual distribution in space of swimming speed could alone explain the observed phenomenon. Therefore, we controlled for these aspects of locomotion across the water tank to make sure no clear gradual pattern emerged (S5 Fig).

Beta oscillations in boundary vector cells

About half of the boundary vector cells exhibited rhythmic neural activity. An example of one such cell is presented in Fig 7A–7C. This cell is a boundary vector cell with a gradually decreasing firing pattern with distance from the bottom of the water tank (Fig 7A). Examining the histogram of interspike interval of this cell (Fig 7B) revealed a periodic spacing pattern between the peaks of the histogram. This was further manifested in the frequency domain: After the histogram was normalized by the total number of spikes, we calculated the power spectral density function of the normalized histogram (Fig 7C), where a local maximum appeared at around 16Hz; in other words, this neuron oscillated rhythmically in the low-beta frequency range. A counter-example of a cell that was not tuned to space is presented in Fig 7D–7F.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Beta oscillations in boundary vector cells.

(A-C) Example of a boundary vector cell with a periodic interspike interval (ISI) pattern. (A) Firing rate heatmap of a boundary vector cell tuned to the bottom of the water tank, color coded from dark blue (zero firing rate) to dark red (maximal firing rate, indicated on the upper right side of the panel). (B) ISI histogram of the cell in A. Even spacing between the peaks suggests periodic oscillations of the neural activity. (C) Power spectral density of the histogram in B (normalized). Colored background marks the typical frequency range of beta waves (12.5–30 Hz). A local peak is shown at approximately 16 Hz, suggesting that this cell exhibited beta-rhythm oscillations. (D-F) A counter-example of a cell with no clear spatial tuning and no specific pattern in the ISI histogram. (G) Maximal power spectral density in the beta waves range for the entire population. The color bar spans the population’s max-p values (see Materials and methods) on a logarithmic scale. As shown, roughly half of the boundary vector cells (blue dots) exhibited beta rhythm oscillations at 15.25 ± 1.62 Hz (mean ± standard deviation). (H) Different thresholds were tested to estimate the prevalence of beta oscillations in the recorded population (see Materials and methods). The median value of the population’s max-p was used to divide the data into groups similar in size. Regardless of the threshold tested, the cells that were more spatially tuned (i.e., below median-max-p = 0.11, blue curve) were more abundant above the threshold than others (i.e., above median-max-p value, orange curve). The underlying data supporting all panels in this figure can be found in a file named Fig 7_data.mat (see Data Availability).

https://doi.org/10.1371/journal.pbio.3001747.g007

Examining the position and magnitude of the peak power spectral density in the beta range of the entire population (Fig 7G, color-coded by the max-p values; see Materials and methods) revealed a clear connection between the beta oscillations and the boundary vector cells in the central telencephalon of the goldfish. To estimate the properties and prevalence of the beta oscillations in the population, we used a threshold of 0.0015 dB/Hz for the peak power of the spectral density. The findings showed that 16 of the 35 boundary vector cells (approximately 46%) and 1 out of the other 84 cells crossed this threshold with a peak spectral density in 15.25 ± 1.62 Hz (mean ± standard deviation). Different thresholds were tested to verify that these results were independent of the chosen threshold (Fig 7H; see Materials and methods).

Sigmoidal firing patterns in space

To further characterize the spatial firing patterns observed in the boundary vector cells, we tested two tuning models: (1) Sigmoid tuning along the cell’s preferred direction, resembling distance encoding; and (2) Gaussian tuning in the form of a 2D Gaussian centered in proximity to the environmental boundaries, resembling local encoding. We fitted the two models for the firing rate map of each cell and tested the goodness of fit for the two models.

This process was first validated using a simulated dataset (examples are presented in S4 Fig; see Materials and methods). Using the simulated data, we set the thresholds (S4J and S4K Fig; see Materials and methods) for classifying a cell as distance encoding or Gaussian encoding (as in a place cell). As expected, the results showed that most boundary vector cells fit better into the sigmoid model than the Gaussian model (see S1 Table). This emerged in the correlations between the data and the models (S4K Fig).

Since a 2D Gaussian with a very broad distribution on only one axis is similar to a sigmoid distribution, these two models are not mutually exclusive. Nevertheless, this analysis helped us to rule out the hypothesis that the cells we recorded were a subset of other typical locally activated space-encoding cells found in vertebrates, characterized by a narrow Gaussian distribution that happen to be centered near boundaries.

Stability of the analysis

To ensure the validity of our analysis over an inhomogeneous swimming trajectory across the water tank (S6 Fig), we ran our analysis pipeline on a set of simulated data. The simulated cells’ activity was drawn from three different spatial tuning curve types (500 shuffled spike trains for each tuning curve). We used a fish’s recorded trajectory, which did not cover the water tank homogenously (corresponding to the trajectory of the cell presented in S2D Fig). Tuning curves of firing rate to position were chosen from three different fish. They were gradually decreasing (S6A Fig, right panel), quasi-uniformly distributed (S6B Fig, right panel), or gradually increasing (S6C Fig, right panel) with the position of the fish along the vertical axis of the water tank (see Materials and methods).

For the decreasing firing pattern, 499 of the 500 simulated cells were indeed classified as spatially modulated. This was also true for 498 out of the 500 cells obtained for the increasing firing pattern (false positive rate <0.004; S6D Fig). None of the cells obtained from the quasi-uniform firing rate (S6D Fig) were classified as spatially modulated (false negative rate <0.002). This confirmed that the analysis for spatial tuning was not sensitive to occupancy correction artifacts in the firing rate map due to partially visited parts of the water tank.

Experimental model

A total of 15 goldfish (Carassius auratus) 13 to 18 cm in length and weighing 100 to 250 g were used in this study. Each fish was housed in a home water tank at room temperature and brought to the experimental water tank for the recording sessions. The room was set to a 12/12-h day–night cycle under artificial light. All the experiments were approved by the Ben-Gurion University of the Negev Institutional Animal Care and Use Committee and were in accordance with government regulations of the State of Israel.

Wireless electrophysiology

The implant preparation and transplanting process, as well as the behavioral electrophysiology details, are fully described in our previous publications [37,38]. Briefly, extracellular recordings were obtained from the goldfish’s central telencephalon using 3 or 4 tetrodes. The data were logged on a small recording device (Mouselog-16, Deuteron Technologies, Jerusalem, Israel) placed in a waterproof case mounted on the fish’s skull. In addition to the tetrodes, a reference electrode was placed near the brain to detect possible motion artifacts rather than neural activity. The tetrodes were moved in the brain between recording sessions using the built-in microdrive of the implant. Wireless communication via a PC and a transceiver (Deuteron Technologies, Israel) was used to control and synchronize the data logger. The neural signals were recorded at 31,250 Hz and band-pass filtered at 300 to 7,000 Hz. A colored styrofoam marker attached to the waterproof case was used to neutralize the buoyancy of the implant (i.e., total average density of 1 g/cm³) and to detect the fish’s position and swimming direction on the video recordings.

Surgery and stereotaxic procedure

During surgery, the anesthetized fish was placed in a holder on the operating table and perfused through its mouth with water and anesthetic (MS-222 200 mg/l, NaHCO3 400 mg/l 1, Cat A-5040, and Cat S-5761, Sigma-Aldrich, USA). All surgery specifics are described in detail in [38]. We targeted the central area of the telencephalon by taking the anterior mid margin of the posterior commissure as the zero point for the stereotaxic procedure [39]. Then, using a mechanical manipulator, we moved the tetrodes 1 mm posteriorly and 1.5 mm ventrally. We used the microdrive attached to the tetrodes to further adjust them in the dorsal/ventral axis after surgery.

Experimental water tank

The experimental water tank was 0.7 m in the horizontal axis and 0.7 m in the vertical axis. The third axis was foreshortened (0.2 m), creating a quasi-2D environment (Fig 1A). The water level during the experiments was constant—0.7 m.

The walls on 3 sides of the tank (one horizontal side and the two short, degenerated sides) were covered and blocked the visual cues from the room. The water tank faced a corner of a room, and the only visual cues from the outside came from the uncovered horizontal side, including a Raspberry Pi and its small camera facing the center of the tank, and a PC (not in front of the tank) communicating with the neural logger (PC screen was covered during the experiments). The top of the water tank was not covered thus ceiling was also visible for the fish.

Video recording

A Raspberry Pi camera (Raspberry Pi 3B microprocessor) oriented towards the center of the water tank was used to localize the fish in the vertical and horizontal dimensions. To synchronize the neural activity and the video recordings, we used Arduino Uno, which sends simultaneous TTL pulses to the data logger’s transceiver and the Raspberry Pi unit. Each recording session lasted 60 to 75 min while the fish navigated freely in the water tank.

Histology

The brains of 11 out of 15 fish were fixed in 4% paraformaldehyde overnight (Electron Microscopy Sciences, CAS #30525-89-4) and then were immersed in a 40% glucose solution for cryoprotection. After freezing, the brains were cryo-sliced (40 μm slices) and Nissl stained to reveal the position of the electrodes in the brain (i.e., the central area of the telencephalon). All slices were then scanned by automated microscope (Panoramic scanner, 3DHISTECH, Hungary) with a 20× (NA 0.8) and a 40× (NA 0.95) objective. The slices were compared to the neuroanatomical landmarks of [40] to confirm the recording region.

Spike sorting

The raw neural data were band-pass filtered at 300 to 7,000 Hz. Then, a threshold was set to detect the action potential timings. Manual single cell clustering was done for each tetrode separately by standard spike sorting methods [41,42], including PCA analysis of the spike amplitudes, widths, and waveforms across all four electrodes. Units not clearly separated in the PCA space and units unstable over time were eliminated from further analysis. In addition, spikes that appeared simultaneously in more than one tetrode or reference electrode were considered motion artifacts or electronic noises and were removed from the analysis. Examples of clean spikes are presented in Fig 1C–1E. More details can be found in [37,38].

Fish trajectory and firing rate-map analysis

The fish’s location and orientation in each video frame were detected using the open-source pose-estimation algorithm by DeepLabCut [43]. This deep-learning algorithm detected several coordinates of the implant on the fish’s head in each video frame. Next, the coordinates were translated to the vertical-horizontal position in cm units. Then, the fish’s trajectory and spike positions were binned using a 5 cm × 5 cm grid and smoothed using a 2D Gaussian filter (σ = 20 cm) to obtain auxiliary maps of occupancy per bin and spike count per bin, respectively. A bin-by-bin division of spike counts by occupancy yielded the occupancy-corrected firing rate map for each cell. Examples of the steps in this correction process are presented in S7 Fig. The corrected rates in each map were color-coded from zero (dark blue) to the maximal firing rate of each cell (dark red). Bins in which the fish spent fewer than 2 seconds were discarded from the analysis and colored white in the figures.

Spatial tuning statistics

To evaluate spatial modulation, we applied the standard statistical tests for spatial modulation in the literature: in-session stability, spatial coherence, and spatial information [44]. For the spatial information, we also used a threshold of 0.1 bits/spike.

To test for in-session stability, we calculated the firing rate maps for the first and second halves of each experiment and calculated the correlation coefficient between them. Then, 5,000 shuffled spike trains were obtained by calculating the interspike intervals (ISIs), shuffling the ISIs by random permutation, and using the cumulative sum to obtain a shuffled spike train. For each shuffled spike train, we created a heatmap using the trajectory of the second half of the experiment and calculated the correlation coefficients between the shuffled heatmap and the heatmap of the first half of the experiment. Comparing the shuffled results with the data yielded a p-value for stability representing the probability of obtaining greater in-session stability than that of the cell by chance.

For the shuffled spike trains, we also calculated the spatial coherence and spatial information values and compared them to the nonshuffled cellular result to obtain two additional p-values. We assigned each cell a max-p value (as used in Figs 6B–6G, 7G and S4K), which was the largest p-value of those calculated in the three tests. Cells with a max-p < 0.025 and for which the spatial information was at least 0.1 bit/spike were counted as spatially modulated cells. This category was used for the purposes of determining the tuning width, rate stability, and directional preferences in spatially modulated neurons and to help estimate the prevalence of beta rhythm oscillations and of boundary vector cells in the population. In the latter case, we also tested different max-p and spatial information thresholds and observed that the prevalence of boundary vector cells was not strongly affected (see S2 Table).

Directional preferences

As one measure of neuronal tuning, we calculated the correlation coefficient between the cellular firing rate and the position in both the vertical and horizontal axes of the water tank. We used a shuffle analysis to get a perspective on the magnitude of these correlations. Specifically, 5,000 shuffled spike trains were obtained by ISI shuffling as described above. For each shuffled spike train, the correlation coefficients were calculated as with the original data. The results of this analysis were used solely for visual assessment (Fig 3A, 3D, 3G and 3J), whereas the actual analysis of the strength of the correlation was taken from the correlation in the preferred and orthogonal directions as defined below.

To determine the preferred boundary direction of the boundary vector cells, we calculated the correlation coefficients between the firing rate and projection of position along 180 axes evenly spaced on the half circle and determined the one for which the value of the correlation coefficient was maximal (Θ_max, the preferred boundary direction) and the one for which the absolute value of the correlation coefficient was minimal (Θ_min). As expected, Θ_min was tightly distributed around the orthogonal direction of Θ_max. To simplify the analysis, we used Θ_max±90^o as the null direction for all neurons.

Then, we calculated the tuning curve of firing rate versus distance from the boundaries found at the Θ_max and Θ_max±90^o directions (as is Fig 3C, 3F, 3I and 3L). This was done by dividing the projection of spike locations on the specified axis by the occupancy using bins of 10 cm. Error bars were calculated by repeating this process 1,000 times with only 50% of the data each time and then calculating the standard deviation in each bin.

Allocentric swimming direction tuning analysis

To test how allocentric swimming direction affected the activity of boundary vector cells, we bisected the data using the displacement direction of the fish. Thus, displacement direction in the range Θ_max±90^o was considered swimming towards Θ_max, and else was considered swimming away it. We used this bisection of the data to plot the heatmaps, tuning curves, and speed distributions shown in Fig 4A–4L.

To quantify the effect of the allocentric swimming direction, we calculated a modulation ratio index (Fig 4M); namely, the ratio that minimizes the least square error between each cell’s two tuning curves—firing rate versus distance from the preferred boundary while swimming towards it and while swimming away from it (e.g., blue and yellow curves in Fig 4B, respectively). For simplicity, this ratio was unsigned, i.e., the stronger tuning curve was divided by the weaker for all cells.

Since behavior might differ while swimming either towards a boundary or away from it, and as an additional measure of stability of this analysis, we calculated the swimming speed distributions of all boundary vector cells along the distance from the preferred boundary for the two bisected data sets (examples are presented in Fig 4C, 4F, 4I and 4L).

Further, we wanted to assess statistically whether the modulation ratio effect might be partially mediated by different swimming speeds towards the boundary versus away from it. Therefore, we used 2 × 2 ANOVAs for each fish with allocentric swimming direction (towards versus away) and distance from the wall (10 cm or 20 cm, i.e., in each cell’s receptive field) as independent variables and speed as the dependent variable. To avoid time correlations in the dependent variable, we down sampled the data to at least 3 seconds between samples: the autocorrelation of the speed goes to 0 at this time difference. This gave us on average 94 independent data points (range: 22 to 322) for the ANOVA for each fish. We determined that a fish swam faster in one direction compared to the other it in case the Bonferroni-corrected p-value of the main effect of direction on speed was at p < 0.05. We determined that the fish swam nearly the same speed towards and away from the boundary if the main effect of direction had a p > 0.15. For other p-values between these two ranges, we classified the effect of direction on speed as indeterminate.

Last, to estimate the relationship between the directional and spatial information in boundary vector cells (Fig 4N), we repeated the process described in [12] to calculate the mutual information between firing rate and location and compared it to the mutual information between firing rate and the allocentric swimming direction of the fish. To do a fair comparison, we used exactly 48 bins (7.5 degrees each) to calculate the directional information and 48 bins on average for the spatial information.

Varying geometry experiment

To test the hypothesis that cells that are tuned to boundaries are also tuned to an additional unfamiliar boundary, in some of the experiments, we altered the geometry of the experimental environment. This was done by placing an additional boundary—a nontransparent horizontal shelf (0.35 m in size), about 0.35 m deep near the right wall of the experimental water tank (see Fig 5B for schematics of the altered environment). This manipulation was done between two subsequent recording sessions while the fish was blindfolded. We recorded 19 of the 119 cells reported in this study in the two different environments.

Beta oscillations analysis

To test for rhythmic oscillations in the beta frequency range, we calculated the histogram of the ISI of each cell using bins of 10 ms. Then, we normalized the histogram and calculated the power spectral density function of the normalized histogram with bins of 1 Hz. For each cell, we calculated the frequency at which the power spectral density was maximal in the beta wave range (12.5 to 30 Hz). Cells with a maximal power spectral density greater than 0.0015 dB/Hz were labelled oscillatory to estimate the prevalence of beta oscillations in the population of boundary vector cells.

To verify that choosing 0.0015 dB/Hz as a threshold did not strongly affect the results, we tested different thresholds (Fig 7H). We used the median-max-p value of the population (0.11) to divide the data into two groups similar in size: those with a max-p value above 0.11 (i.e., with a relatively weak spatial tuning, orange curve) and those with a max-p value below 0.11 (i.e., with a relatively strong spatial tuning, blue curve). The results suggested that regardless of the threshold used, the cells tuned to space (i.e., below the median-max-p value) were more abundant above the different thresholds.

Distance versus local encoding test

Since it remained unclear whether there are “place cells” in the fish telencephalon, and if so, what their properties are, we tested three alternative models to determine what best described the firing patterns observed for the spatially modulated cells: a sigmoid model assuming boundary-vector modulation (i.e., a gradual code along the cell’s preferred direction and a uniform code along the orthogonal direction), a Gaussian model assuming local modulation (as in place cells), and a uniform model of no modulation (control). We tested our ability to differentiate these alternatives on simulated data (S4A–S4J Fig). Each model was tested on a dataset simulated with three different trajectories and three different tuning curves of firing rate to position as recorded in our data.

For the sigmoid model, we simulated 500 gradual firing patterns generated from the three tuning curves of firing rate to position. These were rotated to random preferred directions in space, while in the null (orthogonal) direction, the tuning of rate to position was bounded to a uniform distribution. An example is presented in S4A Fig. Each simulated firing pattern was then fitted with a 2D Gaussian (S4B Fig) and with a sigmoid (S4C Fig).

We used the Matlab fit function to minimize the least square difference from the sigmoid: where L and U are the lower and upper asymptotes of the sigmoid, respectively, x_mid is the point where the curve goes through 50% of the transition, and λ is the width in which 80% of the transition occurs (also used as the boundary vector cells’ width of tuning shown in Fig 6A). The parameter x represents the simulated distance from a boundary in the preferred direction.

The quality of the fit was assessed as the correlation coefficient between the rate map of the simulated data and the rate map generated by the fit model. Each simulated dataset was also fitted with a Gaussian using the Matlab fitgmdist function, and goodness of fit was assessed in the same way.

Simulated firing patterns for the Gaussian model were generated by choosing Gaussians whose centers were within 10 cm of the edges of the environment. This simulated the near border tuning present in our recordings and provided an alternative to the sigmoid tuning. The same pipeline was used for this model as for the former. An example of a simulated Gaussian encoding cell and rate maps fitted from the two models is presented in S4D–S4F Fig. We also repeated the process using 500 rate maps generated from a uniform distribution (S4G–S4I Fig).

The results of the simulation (S4J Fig) showed that the cells exhibiting a gradual firing pattern along a preferred direction could be distinguished from those exhibiting Gaussian firing patterns using the correlation coefficients between the cell’s rate map and the rate maps from the fit models. For this purpose, we drew threshold lines, which created a classification criterion with very few false positives. That is, every neuron with a correlation to the sigmoid model above 0.7 and above the correlation to the Gaussian model was classified as distance encoding. Only 4 of the 500 Gaussian-encoding simulated cells and 7 of the non-encoding (uniform) simulated cells were classified in this way as distance encoding (a false positive rate of around 0.011). In addition, 36 of the 500 cells simulated with a gradual firing pattern were not classified as distance encoding (a false negative rate of around 0.072). We then applied the thresholds derived from the simulations to the real data in order to characterize the firing patterns of the actual neurons (S4K Fig).