Introduction

Foraging in nature entails context-appropriate action selection to resolve a mixture of problems such as risk and value assessment, maintenance of goal-oriented navigation, and task-switching (Pyke et al., 1977; Stephens & Krebs, 1986). For example, foraging animals must update the value of the goal in a space-time domain, negotiate a threat while attaining food, and regulate a variety of stimulus-driven defensive behaviors (Fanselow & Lester, 1988; Kim et al., 2018; Mobbs et al., 2007, 2020; Mobbs & Kim, 2015). All these survival-related problems require executive functions that can flexibly update risk and value representations in working memory to adapt to the foraging environment. Numerous studies have identified the prefrontal cortex as the primary computational circuit for hosting multiple strategies as well as relevant survival-related information (Fuster, 1973; Korn & Bach, 2019; Lee et al., 2014).

A critical component of successful foraging is goal-directed navigation. Although the hippocampus (HIPP) is crucial for maintaining reward location and path planning (Jarzebowski et al., 2022; Krishnan et al., 2022; Nyberg et al., 2022; Wikenheiser et al., 2013; Wikenheiser & Redish, 2015), the medial prefrontal cortex (mPFC) and its interaction with the HIPP also play a role in goal-directed navigation (Guise & Shapiro, 2017; Sapiurka et al., 2016; Spellman et al., 2015; Yu et al., 2018). Furthermore, given that the mPFC is involved in processing deliberate defensive strategies in multiple situations (see Mobbs et al., 2020 for review), naturalistic foraging where risk and reward coexist may serve as a suitable foundation for exploring adaptive information processing in the mPFC. Notably, mPFC houses spatially selective neurons that represent relative positions in a maze (Kaefer et al., 2020; Yang & Mailman, 2018), goal locations (Hok et al., 2005; but see Böhm & Lee, 2020), and precise Cartesian coordinates during food procurement (Mashhoori et al., 2018).

The mPFC also plays a role in processing emotional stimuli, especially those signaling threats in the environment. In line with its intricate connections with the major limbic areas including the amygdala (Gabbott et al., 2005; Hoover & Vertes, 2007; Vertes, 2004), the mPFC is causally linked to learning and expression of fear responses. For example, the infralimbic cortex (IL) of the mPFC is associated with the extinction of conditioned fear response (Park & Choi, 2010; Sierra-Mercado et al., 2006). In contrast, the activation of the prelimbic cortex (PL) is linked to the expression of fear-conditioned responses (Burgos-Robles et al., 2009; Dejean et al., 2016; M. M. Diehl et al., 2018; Kim et al., 2013). In addition, PL neurons show increased activity in response to survival-related stimuli such as food, artificial predators, or both (Kim et al., 2018). These findings suggest that the critical information for risky decision-making might be represented and processed in the mPFC.

Currently, there is limited understanding of how the brain organizes navigation and behavior selections during naturalistic foraging. It is unclear whether spatial representation and emotional processing are separated at the single neuron level or coexist in overlapping populations of the mPFC neurons. In this study, we devised a naturalistic foraging task which allowed self-initiated shuttling between different locations and require avoiding an unpredictable attack by a robotic predator while approaching sucrose. This approach allowed continuous monitoring of the neural activity as the rat switched between navigation and strategic reward procurement. Machine learning-based populational decoding methods, alongside single-cell analyses, were employed to investigate the correlations between neuronal activity and a range of behavioral indices across different sections within the foraging arena.

Results

Mix of avoidance and escape behaviors during naturalistic foraging under threat

The foraging arena, modified from (Kimm & Choi, 2018), was designed to encourage naturalistic foraging involving both approach to food and avoidance of threat. The arena was divided into two distinctive zones by a wall with a pass-through opening in the middle (Figure 1A): the foraging (F) and nest (N) zones. At one end of the F-zone, a motorized gate controlled the access to the robot compartment, which houses a sucrose port and the predator-like Lobsterbot, named after a pair of its large, claw-like arms. The rat’s entrance to the sucrose port was detected by the beam sensor mounted near the gate (Figure 1B). A semicircular area surrounding the robot compartment was designated as encounter (E) zone where the rat’s body was situated during sucrose licking and Lobsterbot attack.

Behavioral procedure and apparatus.

(A) Schematic illustration of the arena showing major sections (N: nest zone, F: foraging zone, E: encounter zone). At the end of the E-zone, Lobsterbot (red) is situated, guarding the sucrose delivery port (green). (B) A rendered 3-D image of Lobsterbot. The sucrose port is located between the “claws”. The two red lines indicate infrared detectors, one for lick detection (short line) and the other for the entry to the E-zone (long line). (C) The experimental schedule. (D) Sample snapshots of Avoidance Withdrawal (AW) and Escape Withdrawal (EW). In an AW trial, the rat typically retracts its head ahead of time and watches the Lobsterbot attack. In an EW trial, the rat reflexively flees from the attack. (E) Example behavior data containing two consecutive trials (Trial 15 and 16). Each trial started with a reentry to the N zone which triggers gate opening. The rat leaving the N zone typically moves toward the E-zone across the F-zone. The entry to the E zone is detected by an IR beam sensor (blue shade). Within the E zone, the rat starts licking (green lines) until being attacked by Lobsterbot (red line) 3 or 6 s after the first lick. The rat shows voluntary withdrawal behavior (AW; Trial 15) or forced escape behavior (EW; Trial 16). (F) A summary of the AW trial rates for each animal during the Losterbot sessions. Points for Lob2 of Rat2 and Rat3 are omitted because they did not approach to robot during the entire Lob2 session.

The experiment was composed of three distinct phases (Figure 1C): habituation phase, shuttling phase, and Lobsterbot phase. The Lobsterbot phase was further divided into early and late phases separated by the time of surgery. During the habituation phase, a food-deprived rat was acclimated to the sucrose port that released sucrose solution. The lick port was activated by an IR-beam sensor, triggering the solenoid valve when the beam was interrupted. The rat gradually learned to obtain rewards by continuously licking the port. Once the rat maintained a stable licking frequency, typically after 3 days, the next phase started. During the shuttling phase, the rat was trained to shuttle back and forth to procure sucrose. The robot compartment’s gate closed 6 seconds after the first lick and reopened only when the rat returned to the N-zone. The shuttling phase lasted 3 days, followed by the introduction of the Lobsterbot phase. During this phase, the robot executed a fast-striking attack, 3 s (30%) or 6 s (70%) from the first lick, with the timing randomized across trials to mimic natural uncertainty and prevent simple time-counting strategies. As in the shuttling phase, the gate was automatically closed after the attack, and the rat must return to the N-zone to reopen the gate. In each trial of the Lobsterbot phase, the foraging rat could respond with either an avoidance withdrawal (AW: retraction of the head before the attack) or an escape withdrawal (EW: reflexive retraction after the attack; Figure 1D). Figure 1E shows a representative sequence of responses and events in two consecutive trials. The AW was clearly distinguished by a wide gap between head withdrawal (indicated by the end of the blue shade) and the attack (red line). Conversely, the characteristic feature of the EW was the overlap between the attack and the animal’s presence in the E-zone (Video 1 & 2).

We monitored all head withdrawal timepoints to assess whether rats developed a temporal strategy to differentiate between the 3-second and 6-second attacks. We found no evidence of such a strategy, as the timings of premature head withdrawals during the 6-second attack trials were evenly distributed (see Supplementary Figure S1).

All five rats exhibited distinctive phasic behavioral changes. Notably, there was a drop in approach attempts (i.e., number of trials) during the first session of the Lobsterbot phase (Lob1), followed by a stable rate of approach attempts in subsequent sessions after the surgery (starting from Lob4) (Supplementary Figure S2A). Other response measures, such as the number of licks (Supplementary Figure S2B), licks per trial (Supplementary Figure S2C) and lick latency after gate opening (Supplementary Figure S2D), also showed a similar phasic change. To assess these behavioral changes, we compared the last session of the shuttling phase (Sht3) and all sessions of the Lobsterbot phase before the surgery (Lob1, Lob2, Lob3). The analysis showed that the number of approach attempts (Friedman test; χ²(3)=9.240, p=.017), the number of licks (χ²(3)=12.60, p <.001), the lick latency (χ²(3)=7.800, p =.044), and the licks per trial (χ²(3)=13.08, p <.001) were significantly changed during these sessions. More specifically, number of licks and licks per trial were decreased at the first session of the Lobsterbot phase (Dunn’s multiple comparisons test; p=.021; p=.043) and subsequently increased back at the end of the early Lobsterbot phase (Lob3; p>.999; p>.999). In summary, after the three Lobsterbot sessions, the rats’ number of approach attempts, number of licks, licks per trial, and lick latency reverted to shuttling phase levels.

In the early sessions of the Lobsterbot phase, the rats showed passive defensive behaviors such as freezing (Blanchard & Blanchard, 1969; McClelland & Colman, 1967) and stretched attend posture (Grewal et al., 1997; Mackintosh & Grant, 1963). However, these behaviors gradually diminished toward the middle of the Lobsterbot phase. Compared to the shuttling phase, during which the rats continue to lick the sucrose port until the gate is closed on almost all trials (88.96%), a mix of AW and EW (65.11% and 34.89% respectively) trials were observed during the Lobsterbot phase. Due to the transient nature of the decreased approach attempts and the various defensive behaviors after the first encounter, surgery was performed after three Lobsterbot sessions (Figure 1F).

Population activity in the mPFC predicts the distance from the Lobsterbot while navigating in the arena

Recording data were collected from a total of five rats implanted with a custom-made moveable tetrode drive. The tetrodes were initially implanted in the PL and advanced ventrally toward the IL at the end of each session (0.1∼0.2 mm/session). Recording locations were reversely calculated from the marking lesion made after the last recording session. Tetrode tracks were identified by visual inspection under a microscope as indicated by black vertical lines on the matching atlas (Paxinos & Watson, 2007) as shown in Figure 2A. All tetrodes except for one spanned both the PL and the IL. After cluster cutting and preprocessing (see Materials and Methods), 632 units (57 ∼ 198 units per rat) from the PL (n = 379) and the IL (n = 253) were selected for ensemble analysis.

Ensemble activity from the mPFC predicts distance from the goal.

(A) Schematic diagram of electrode implantation and estimated recording site. Top: A movable 4-tetrode microdrive was initially implanted in the PL region and lowered ventrally toward the IL after every recording session. Bottom: Representative recording tracks from all five animal are superimposed over an image of a stained coronal section of the frontal brain. Histological examination of all brain sections confirmed that the electrode tracks spanned the dorsoventral axis between the PL and IL. (B) Modulation of unit firing showing place-cell like activities. Units 66 and 125 exhibit fragmented place fields in all over the arena, while Units 56 and 26 display relatively large place fields surrounding particular spots such as the gates. Heat maps are calculated from z-scored spatial tuning curves. (C) Schematics of the ensemble decoding analysis. The 4-layer deep artificial neural network (ANN) receives populational neural data during 50ms-timewindow and is trained to predict the rat’s current distance from the center of the E-zone. The example data depicted in the figure is a sample recording from 20 units when the rat is at a particular distance away from the center of the E-zone, indicated by the white bold line. (D) Accuracy of the regressor. Mean Absolute Error (MAE) was calculated for the two types of regressor: one with the original dataset (Original) and another with the shuffled dataset (Shuffled). The average MAE was 16.61 cm for the Original, which was significantly smaller than that for the Shuffled and smaller than the rat’s body size. This suggests that the mPFC might encode the spatial correlates reflecting distance from the goal. (E) Prediction accuracy in the F-zone during outbound/inbound paths. Decoding accuracy in the F-zone was calculated separately for the outbound (from the N-zone to the E-zone) and inbound (from the E-zone to the N-zone) paths. The decoding accuracy remained unchanged despite the differences in motivation and perceived visual cues due to the movement direction. (F) Comparison of the regressor’s accuracy from the control experiment. When the Lobsterbot was removed from the robot compartment, reverting the task back to simple shuttling, the mPFC distance regressor’s performance significantly decreased compared to the Lobsterbot phase.

Previous studies have reported spatial correlates were observed in multiple subregions of the mPFC (Hok et al., 2005; Jung et al., 1998; Ma et al., 2023; Mashhoori et al., 2018; Sauer et al., 2022; Spellman et al., 2015). To determine whether location-specific neural activity exists at the single-cell level in our mPFC data, a traditional place map was constructed for individual neurons. Although most neurons did not show cohesive place fields across different regions in the arena, a few neurons modulated their firing rates based on the rat’s current location. However, even these neurons were not comparable to place cells in the hippocampus (O’Keefe & Dostrovsky, 1971) or grid cells in the entorhinal cortex (Hafting et al., 2005) as the place fields were patchy and irregular in some cases (Figure 2B; Units 66 and 125) or too large, spanning the entire zone rather than a discrete location within it (Units 26 and 56). The latter type of neuron has been identified in other studies (e.g., Kaefer et al., 2020).

Since we did not find neural activity that can be decoded into accurate spatial information based on individual activity, we performed machine-learning (ML) analysis targeting all recorded units during the given session. Among different ML algorithms, we implemented a four-layer artificial neural network regressor (ANN; see Materials and Methods for a detailed structure), as the ANN achieves significantly lower decoding errors (Supplementary Figure S3) and has robustness to hyperparameter changes (Glaser et al., 2020). The inputs to the regressor were the simultaneously recorded neural activity which was segmented into 50-ms bins and the predicted output was the distance between the rat’s head and the center of the E-zone where reward and threat co-exist (Figure 2C).

The analysis of trained ANNs revealed that the distance from the rat to the center of the E-zone can be decoded from the mPFC ensemble activity. Compared to the control regressor trained on a shuffled dataset, the regressor with original data showed significantly lower errors (paired t-test, t(39) = 10.47, p <.001). However, we did not observe a difference in decoding accuracy between the PL and IL ensembles, and there were no significant interactions between regressor type (shuffled vs. original) and regions (mixed-effects model; regions: p=.996; interaction: p=.782). These results indicate that the population activity in both the PL and IL contains spatial information (Figure 2D, Video 3). The overall accuracy defined in the form of the mean absolute error (MAE) was 16.61 (± 3.67) cm, which was slightly higher than values from previous reports employing different spatial tasks (Ma et al., 2023; Mashhoori et al., 2018). In the current study, an average of 15.8 neurons, a number lower than that of previous studies, was included in a given ensemble in our analysis. Considering the nature of the machine-learning algorithm, we expected that accuracy would be higher with a greater number of simultaneously recorded neurons.

To determine whether the accuracy of the regressor varied depending on the direction of movement, we compared the decoding accuracy of the regressor for outbound (from the N-to E-zone) vs. inbound (from the E- to N-zone) navigation within the F-zone. There was no significant difference in decoding accuracy between outbound vs. inbound trips (paired t-test; t(39) = 1.52, p =.136), indicating that the stability of spatial encoding was maintained regardless of the moving direction or perceived context (Figure 2E).

Previous analyses indicated that the distance regressor performed robustly regardless of movement direction, but there is a possibility that the decoder detects visual cues or behaviors specific to the E-zone. For example, neural activity related to Lobsterbot confrontation or licking behavior might be used by the regressor to decode distance. To rule out this possibility, we analyzed a subset of data collected when the compartment door was closed, preventing visual access to the Lobsterbot and sucrose port and limiting active foraging behavior. The regressor trained on this subset still decoded distance with a MAE of 12.14 (± 3.046) cm (paired t-test; t(39) = 12.17, p <.001). Notably, the regressor’s performance was significantly higher with this subset than with the full dataset (paired t-test; t(39) = 9.895, p <.001). In summary, these data suggest that precise spatial information could be extracted using the neural ensemble in the mPFC.

Removing the Lobsterbot decreases spatial decoding accuracy

Earlier studies have suggested that PL neurons do not exhibit any spatially correlated firing patterns in an empty round apparatus (Poucet, 1997). This absence of spatial correlates may result from a lack of complex goal-oriented navigation behavior, which requires deliberate planning to acquire more rewards and avoid potential threats. To evaluate how the Lobsterbot’s presence impacted spatial decoding accuracy, additional single-unit recording experiments (Ctrl-Exp) were conducted and compared with the previous results (Lob-Exp). In this experiment, three rats followed the same protocol until the microdrive implant surgery. After the surgery, unlike the Lob-Exp group, the Ctrl-Exp group returned to the shuttling phase, during which the Lobsterbot was removed. With this protocol, both groups experienced sessions with the Lobsterbot, but the Ctrl-Exp group’s task became less complex, as it was reduced to mere reward collection. During a total of nine recording sessions, 112 units were identified and ANN regressors were built and trained using these data. Statistical tests indicated that the original dataset had significantly lower decoding errors (20.48 ±2.31 cm) compared to the shuffled dataset (22.61 ±1.95 cm; paired t-test; t(8) = 6.20, p <.001). Nevertheless, when the error was compared with that of Lob-Exp, the Ctrl-Exp group showed a higher decoding errors (Figure 2F; unpaired t-test; t(47) = 3.02, p =.004). This did not result from differences in activity levels, as there was no significant difference in total travel distance between the Lob-Exp and Ctrl-Exp groups (unpaired t-test; t(47) = 0.07, p =.941).

Another factor potentially contributing to the increase in decoding errors could be differences in recording quality across experiments. To discriminate the effects, a generalized linear model (GLM) was constructed with predictors of 1) the Lobsterbot’s presence, 2) the number of units recorded, and 3) the recording site (PL vs. IL) on decoding error. The results showed that the GLM predicted the outcome variable with significant accuracy (F(3,45) = 14.06, p <.001), explaining 48.38% of the variance in decoding error. The regression coefficients were as follows: presence of the Lobsterbot (2.76, standard error [SE] = 1.11, t = 2.49, p =.016), number of recorded cells (−0.43, SE =.08, t = 5.29, p <.001), and recording location (0.91, SE =.92, p =.329). The results indicated that the number of recorded cells significantly influenced decoding accuracy. Most importantly, the presence of the Lobsterbot significantly decreased the overall decoding error, even after accounting for the number of recorded cells. In summary, the removal of the Lobsterbot decreased spatial decoding accuracy.

Non-navigational behavior reduces the accuracy of decoded location

Although initial analysis showed that the distance of the rat to the center of the E-zone could be decoded with reasonable accuracy, the spatial precision was unevenly distributed throughout the arena (Figure 3A). The magnitude of error was highest in the N-zone, followed by the E-zone, and lowest in the F-zone. A one-way repeated measures ANOVA revealed a significant main effect of the zone (F(1.68, 65.52) = 40.36, p <.001). Post hoc analysis showed that the MAE in the F-zone (13.86 ± 2.56 cm) was significantly lower than that in the N-zone (22.72 ± 6.38 cm; t(39) = 10.69; p <.001) and the E-zone (19.98 ± 7.30 cm; t(39) = 6.34; p <.001). However, there was no significant difference between the N- and E-zones (t(39) = 2.29; p =.081; Figure 3B). To correct for the unequal distribution of location visits (more visits to the F-than to other zones), the regressor was trained using a subset of the original data, which was equalized for the data size per distance range (see Materials and Methods). Despite the correction, there was a significant main effect of the zone (F(1.16, 45.43) = 119.2, p <.001) and the post hoc results showed that the MAEs in the N-zone (19.52 ± 4.46 cm; t(39) = 10.45; p <.001) and the E-zone (26.13 ± 7.57 cm; t(39) = 11.40; p <.001) had a significantly higher errors when compared to the F-zone (14.10 ± 1.64 cm).

Spatial encoding is disrupted by non-navigational behaviors

(A) Spatial distribution of the prediction accuracy. The heatmap indicates MAE of a fully trained ANN, superimposed over the entire foraging arena. The prediction accuracies were lower in the N- and E-zone than that in the F-zone. (B) Mean prediction accuracy by the zones. The MAE in F-zone was significantly lower than the other zones. Error bars represent the SEM. (C) Examples of the non-navigational behaviors in the N-zone. The top three snapshots depict grooming, rearing and sniffing. The bottom three snapshots show typical goal-directed navigational movements. (D) Comparison of decoding errors (N-zone) during navigational vs. non-navigational behaviors. The error was significantly larger when the rat was engaged in non-navigational behaviors within the N-zone. (E) Comparison between regressors trained with and without non-navigational behaviors. The overall decoding error was significantly smaller when the regressor was trained without the data during non-navigational behaviors.

In the N-zone, the rats engaged in stereotypic behaviors such as grooming, climbing, and rearing, during which they stopped locomotive behavior (Figure 3C). To determine whether these non-navigational behaviors resulted in higher decoding errors, these behaviors inside the N-zone is manually labelled and the MAEs during navigational vs. non-navigational behaviors were compared. Due to the imbalanced nature of the dataset resulting from varying preference for non-navigational behaviors, data from 21 sessions were included, each containing at least 20% of one type of behavior. The paired t-test showed a significant difference between the MAEs during the non-navigational and navigational behaviors in the N-zone (paired t-test; t(20) = 9.381, p <.001; Figure 3D). Moreover, excluding the non-navigational data of the N-zone from the original dataset improved the overall accuracy of the regressor (paired t-test; t(20) = 8.562, p <.001; Figure 3E). These data indicate that erroneous decoding result in the N-zone might have resulted from the disruption induced by non-navigational behaviors.

Population analysis reveals a different functional mode in the E-zone

To assess the functional coherence of population activity (Durstewitz et al., 2010) within and across different zones, principal component analysis (PCA) was performed on the ensemble data collected throughout the session. Figure 4A illustrates the first two principal components of the transformed neural data from a representative recording session (number of units = 19). First, within- and between-zone distances in population neural activity were compared to determine whether activity within the same zone exhibited similar properties, distinct from those of other zones A paired t-test revealed a significant difference between distances (t(119) = 10.46, p <.001), confirming that the neural activity within a given zone was distinct from those of other zones.

PCA results reveals distinctive population activity in the E-zone

(A) Representative recording session depicting the first two dimensions of the PCA result. Populational neural activities are projected onto a virtual space. Each dot represents 50 ms-long neural activity from multiple units. The color of the dots indicates the rat’s location during the corresponding neural activity. Diamonds represent the centroids of neural representations for each zone. To visually emphasize each cluster, data points close to centroids are selectively plotted. (B) Distances between each centroid pairs from all recording sessions. The centroid of the E-zone is distinctly positioned compared to the centroids of the other two zones, indicating a unique neural state within the E-zone. The triangle above the bar plot represents the relative distance between the centroids of each zone’s neural ensemble activity. Longer edges signify greater dissimilarity between the neural ensemble activities of two zones. Error bars represent the SEM.

To quantitatively measure functional resemblance among ensemble activity within different zones, we measured the distance between pairs of centroids projected onto the latent space (N-F, N-E, and F-E). A Friedman test showed a significant difference between centroid pairs (χ²(2)=51.05, p <.001). Further post hoc analysis revealed a significant difference between N-F and N-E pairs (p <.001) and N-E and F-E pairs (p <.001) but not in N-F and F-E pairs (p =.353), suggesting an ordinal similarity with N-F (D =.72) being most similar, followed by F-E (D = .94) and N-E (D = 1.44) (Figure 4 Figure 4B).

Not surprisingly, the firing pattern in the E-zone was the most dissimilar, reflecting a functional state different from those in other zones. In contrast, the distance between centroids of the N- and F-zones was relatively small, indicating that the overlap between the neural subspaces in the two zones was not negligible.

Although PCA can capture population differences between neural ensembles in the zones, a small subset of neurons highly tuned to salient sensory cues could solely drive this separation. Particularly in the E-zone, which encompasses two significant events—head-entry and head-withdrawal—neurons responsive to the delivery of rewards or attacks from the robot could dominate the principal components. To minimize this influence, we defined “critical event times” as ±1 second from head-entry and head-withdrawal, excluded neural data during these periods, and then recalculated the principal components. The results still demonstrated the same sequential ensemble differences (F(1.57, 61.34) = 73.56, p <.001), indicating that the population differences are not driven by events specific to any zone (Supplementary Figure S4). In summary, populational activity inside the E-zone showed distinctive activity patterns unlike those of active navigation.

Hierarchical clustering identifies functionally distinctive sub-populations in the E-zone

In the E-zone, navigational behavior was mostly replaced by approach attempts and avoidance behaviors. To characterize the temporal dynamics of single-unit activity in the E-zone, we examined the modulation of firing rate to the two most critical behaviors; the head-entry and the head-withdrawal. The temporal dynamic of neuronal firing relative to these two behavioral events was visualized using color-coded raster plots and peri-event time histograms (Figure 5A & 5C), featuring spiking activities from all behavior-responsive units (see Materials and Methods). Overall, the units showed modulated firing with varying degrees of temporal relations to the two events. Most strikingly, overall modulation patterns were clearly different between head-entry and head-withdrawal.

Multiple subpopulations in the mPFC react differently to head entry and head withdrawal.

(A) Top: The PETH of head entry-responsive units is color-coded based on the Z-score of activity. Bottom: The red vertical lines mark the timing of the head-entry. The peak latency of each unit varies from as early as 2 s before and to 1∼2 s after the head-entry. (B) Functional segregation of all recorded units. Top and middle: Two sub-populations of units based on hierarchical cell clustering analysis. Bottom: The averaged activity for each sub-population. (C) The PETH of head withdrawal-responsive units is color-coded based on the Z-score of activity. (D) Functional segregation of all recorded units. Top and middle: Three sub-populations of units based on hierarchical cell clustering analysis. Bottom: The averaged activity for each sub-population.

For the head-entry, a large portion of the units increased the firing before the event (Figure 5A, top). The pattern was clearly present when all neuronal spike activities were collectively analyzed (Figure 5A, bottom). The mean activity showed a peak approximately 150 ms before the head-entry, with the peak exhibiting substantial width (FWHM of 200 ms), suggesting the presence of varied neuron populations differing in peak position. In contrast, this anticipatory modulation was not clearly identified when the activities were aligned to the head-withdrawal. Instead, there was a strong temporal aggregation of activity around the peri-event period (Figure 5C, top) within the narrow range. This was evident in the collapsed histogram as a single sharp peak at the time of the head-withdrawal (peak at 0 ms, FWHM of 75 ms; Figure 5C, bottom).

To categorize recorded units into multiple sub-populations based on their functional connectivity in reference to the key events, the unsupervised hierarchical clustering method was applied to all recorded units (Boehlen et al., 2016; Di Miceli et al., 2020). We identified two sub-populations when we classified neurons according to their activity around the head-entry (Figure 5B). The first sub-population (HE1: n = 366, 57.91%) showed an anticipatory increase with a peak located at 350 ms before and a negative modulation after entry and remained low for an extended period. The other sub-population (HE2: n = 189, 29.91%) showed a more concentrated positive modulation at around the time of entry. When we used the head withdrawal as criterion, three sub-populations were identified (Figure 5D). The first (HW1: n = 380, 60.13%) showed a complex modulation pattern: lower than the baseline in the beginning, followed by a peak firing timed at the head-withdrawal. The second (HW2: n = 147, 23.26%) maintained increased activity until the onset of head-withdrawal and returned to zero afterward. The last sub-population, HW3, which consisted of a relatively small number of units (n = 62, 9.81%), showed an inhibitory response at the time of head-withdrawal.

Well-trained rats typically run across the F-zone at high velocity and arrive at the sucrose port in the E-zone, stopping abruptly. Moreover, the head-entry is always aligned with the rats’ movements, including both subtle (AW) and reflexive (EW) movement. Consequently, it is plausible that the sudden movement might result sudden change in neural activity around the time of head-entry and/or head-withdrawal. To exclude this possibility, a “run-and-stop” event was correlated with the neural data. The event was defined as the transitional moment between a directional locomotion of 2 s or longer and a subsequent pause of 1 s or longer (see Materials and Methods). Unit activities were aligned around this run- and-stop event, and the mean activity of each sub-population was plotted for visual inspection. As shown in the Supplementary Figure S5, there was no detectable neural modulation around the run-and-stop events. The lack of significant modulation at the time of action suggests that the neural modulation occurring around the head-entry and head-withdrawal is not a simple motoric response.

An interesting overlap was found between sub-populations characterized by head-entry and head-withdrawal. In the HW1 group, 78.68% (n = 299) of units were classified as HE1. On the other hand, 63.95% (n = 94) of the HW2 units were from HE2 (Supplementary Figure S6A). Based on this, we further characterized the temporal dynamics of spike activity in relation to the two behavioral events using the combined categories: Type 1 (HE1-HW1) and Type 2 (HE2-HW2). These two subpopulations demonstrated distinctive firing patterns within the E-zone (Supplementary Figure S6B & S7). Specifically, Type 1 neurons exhibited a rapid decline in activity following a short peak 200 ms before the head-entry, whereas Type 2 neurons displayed a sustained high peak activity (z ≍ 2) from the onset of the head-entry until the head-withdrawal.

We also examined whether there is a significant difference between the PL and IL in the proportion of Type 1 and Type 2 neurons. In the PL, among 379 recorded units, 143 units (37.73%) were labeled as Type 1, and 75 units (19.79%) were labeled as Type 2. In contrast, in the IL, 156 units (61.66%) and 19 units (7.51%) of 253 recorded units were labeled as Type 1 and Type 2, respectively. A Chi-square analysis revealed that the PL contains a significantly higher proportion of Type 2 neurons (χ²(1, 632) = 34.85, p < .001), while the IL contains a significantly higher proportion of Type 1 neurons compared to the other region (χ²(1, 632) = 18.07, p < .001). Taken together, these findings provide support for the existence of at least two distinct subpopulations in the PL and IL that may be involved in different cognitive functions within the E-zone.

Population activities in the E-zone encode success and failure of avoidance with high accuracy

The previous analysis confirmed that a significant number of neurons modulated their activity before or around the time of approach toward and withdrawal from the E-zone. We then further analyzed the correlation between neural activity and success and failure of avoidance behavior. If the mPFC neurons participate in the avoidance decisions, avoidance withdrawal (AW; withdrawal before the attack) and escape withdrawal (EW; withdrawal after the attack) may be distinguishable from decoded population activity even prior to motor execution. A naïve Bayesian classifier (see Materials and Methods) was trained to distinguish AW from EW using 2-s neural ensemble data around the onset of the event (Figure 6A). We corrected the accuracy against the sample imbalance owing to varying degrees of the AW/EW response ratio across different sessions. Our analysis showed that the classifier based on the original dataset achieved a significantly higher prediction accuracy (72.71%) than a theoretical chance-level prediction (50%) or a shuffled dataset (51.62%) (paired t-test; t(39) = 10.06, p <.001) in distinguishing AW from EW (Figure 6B). Furthermore, we analyzed whether there is a difference in prediction accuracy between sessions with different recorded regions, the PL and the IL. A repeated two-way ANOVA revealed no significant difference between recorded regions, nor any interaction (regions: F(1, 38) = 0.1828, p = 0.671; interaction: F(1, 38) = 0.1614, p = 0.690).

Neural ensemble activity predicts failure and success of avoidance response.

(A) Schematics of the event decoding analysis. The Naïve Bayesian decoders are trained with 2 s window of neural activity to discriminate avoidance or escape on every trial (AW/EW classifier). The grayscale image depicts an example firing pattern of 17 units on a given trial, arranged to the onset of the withdrawal response. The decoder classifies whether this trial is AW or EW based on this data. (B) Accuracy of the naïve Bayesian classifier. The decoding accuracy of the classifier was significantly higher than that from the shuffled data. (C) Temporal characteristics of prediction accuracy of the naïve Bayesian classifier. Prediction accuracy was significantly higher at the time points as early as 5-7 s before the head-withdrawal. (D) Class discrimination index by the two sub populations of neurons. The class discrimination index indicates that the Type 2 neurons showed a significant discriminatory power towards AW. Neurons in the Type 2 and the Others group did not exhibit significant discriminatory power.

Because the behavioral reactions following the two events were quite different (bigger reactions during the EW as if the animal failed to predict the attack), it is not clear whether the naïve Bayesian classifier discriminated the two types of responses based on the post-event population activity, which would have been a reflection of the different reactions in the two responses rather than an anticipatory or value-driven encoding. To clarify this issue, we trained naïve Bayesian classifiers using 2-s ensemble activity at varying time points before the execution of the head-withdrawal response. Multiple t-tests between original and shuffled datasets revealed that AW/EW discrimination accuracy began to increase above the chance level 6 s before the actual response and remain significant at all time points (multiple paired t-tests with Šídák-Boferroni correction; ps < .01; Figure 6C, Supplementary Data 1).

Based on the previous clustering analysis that identified functionally cohesive sub-populations of neurons (Types 1 and 2), we further investigated the predictive encoding by the ensemble activity of the two subpopulations. The naïve Bayesian classifier uses class likelihood, P(unit|AW) and P(unit|EW), to classify a particular type of withdrawal response. By comparing the class likelihoods, we can assess the ability of these subgroups to distinguish between AW and EW. This ability is quantified by a ’Class Discrimination Index,’ calculated as a log odds ratio. This index reflects the likelihood that a given feature set is associated with a particular type of withdrawal response (AW or EW). To test whether the unit activity right before the head-withdrawal (from 2 s before the head-withdrawal to the onset of the head-withdrawal) was correlated with subsequent avoidance behavior, we calculated the log odds ratio from the trained classifier. The results showed that Type 2 neurons had a significantly higher log odds ratio, indicating a high correlation between unit activity and success in avoidance, whereas no discriminatory modulation was observed for Type 1 neurons or others (one sample t-test, t(298) = .564; p =.573, t(94) = 5.63; p <.001, t(238) = .531; p =.596, Type 1, Type 2, and others, respectively). These results indicate that Type 2 units carries information that is relevant for classifying trial as AW. To summarize, the populational neural activity of the mPFC in rats within the E-zone distinguished between two types of head-withdrawal behavior, and increased activity of unit subgroups was strongly associated with the execution of AW.

Overlapping populations of mPFC neurons adaptively encode spatial information and defensive decision

Through a series of analyses, we found that both distance and foraging decision can be decoded from the same mPFC neural ensemble. PCA results support the that the mPFC ensemble activity forms a distinctive functional subspace when the rat is engaged in foraging decision facing the robotic predator in the E-zone. Previous studies have also observed different firing patterns from the same mPFC neural ensembles when the rat switches between tasks (Durstewitz et al., 2010; Nigro et al., 2023). Therefore, it would be interesting to see whether the distance and foraging decision are encoded by two exclusively subsets of neurons or by the whole ensemble.

To quantitatively assess each neuron’s contribution to the distance regressor and event classifier, we calculated the feature importance scores for each decoder (Figure 7A). This metric measures the increase in error when a neuron’s data is shuffled. For example, if a neuron holds substantial data and is heavily relied upon by the decoder, shuffling its signal will significantly increase the error. We analyzed the permutation feature importance scores for each neuron for both distance regression and event classification (see Materials and Methods).

Feature importance analysis found no evidence of dedicated neural subset for Distance/Event encoding.

(A) Schematic diagram showing computational protocols for the feature importance analysis

(B) Frequency distribution of all recorded units for their feature importance (measured in error increase) in distance encoding. Only a few units produced non-negligible increase in error when removed from the data indicating the absence of dedicated distance-encoding neurons. Red line indicates 95th percentile.

(C) Accuracy of distance regressor without top 20% of high-performance units. Even when 20% of units with high feature importance score were removed, the distance regressor could decode rat’s location.

(D) Frequency distribution of all recorded units for their feature importance (measured in accuracy drop) in event classification. Only a few units produced non-negligible decrease in accuracy when removed from the data indicating the absence of dedicated event-encoding neurons. Shuffling unit’s data resulted small decrease in error, and even some showed increase in accuracy. Red line indicates 95th percentile.

(E) Accuracy of event classifier without top 20% of high-performance units. Even when 20% of units with high feature importance score were removed, the event classifier could decode the type of rat’s defensive behavior.

(F) Correlation between distance regressor’s feature importance score and event classifier’s feature importance score. There was no correlation between feature importance of two types of decoding

First, for the distance regressor, the median permutation feature importance score was 0.73 cm (Figure 7B). This suggests that shuffling one unit’s neural activity increases the decoding error by less than 1 cm. Moreover, 95% of units exhibited a permutation feature importance score below 1.998 cm, indicating only a negligible increase. To determine whether a subset of neurons with high feature importance scores solely drives distance decoding, we tested whether removing the top 20% of neurons with the highest scores in a given session affects distance decoding. The results showed that although removal significantly increased the decoding error, the regressor trained with the remaining data still achieved a significantly lower error (17.35cm; Wilcoxon signed rank test, p<.001). Next, we applied a similar method to the event classifier. This analysis revealed a median permutation feature importance score of 0.01 and a 95th percentile of 0.06. We also tested the effect of removing the top 20% of neurons with the highest scores and found that the decoder still exhibited a significantly lower error (paired t-test, t(39)=6.899, p<.001). This result indicates that there is no sub-population which dominates in encoding distance information.

Furthermore, we tested whether there is a correlation between the feature importance scores of the distance regressor and the event classifier. If there were two subgroups of neurons, each dedicated to encoding distance and events respectively, we might expect to see a negative correlation. To verify this, we calculated Pearson’s correlation for each session. We found that there was no significant correlation between the feature importance scores of the two encoders (r=0.077, p=.062), indicating the absence of two distinct subpopulations that exclusively encode either distance or event.

Additionally, we found no evidence that Type I or Type II neurons have statistically higher or lower feature importance in the distance regressor. These results support the notion that distance and avoidance outcomes are not encoded separately by subsets of neurons; rather, the entire ensemble encodes two different concepts according to the given context.

Discussion

We found distinctively different states of ensemble activity from the mPFC (PL and IL) of a foraging rat, and these states were linked to specific areas or tasks within the foraging environment. Specifically, the rat’s distance from the conflict zone could be decoded with stable accuracy during navigation, while its ensemble state switched to encode a defensive behavior when the rat entered the immediate vicinity of a reward and threat region (E-zone). Within the high-conflict area, the rat’s threat-evading behavior (AW), coupled with forgoing the reward, could be decoded using a naïve Bayesian classifier. These results demonstrated that the rat’s mPFC neurons were closely associated with two different functional modes of processing allocentric spatial information as well as egocentric information, the reward and threat.

In the present study, the rat’s spatial information was extracted from the mPFC neuron in the form of distance. Previous studies examining the mPFC’s spatial representation found a spectrum of outcomes, from an absence of spatial information (Poucet, 1997) to accurate (x, y) coordinates with errors as small as a few centimeters (Ma et al., 2023; Mashhoori et al., 2018). In some studies, the spatial coding mechanism was directly compared between the hippocampus and the mPFC (Kaefer et al., 2020; Powell & Redish, 2014), and the mPFC cells were found to encode a generalized form of space rather than specific locations. Therefore, we focused on distance as a more generalized spatial index that also reflects an internal representation of value and risk, and built an ANN regressor to target it.

The decoding analysis showed that brief navigation was sufficient to engage a majority of mPFC neurons recorded during foraging into spatial encoding. Furthermore, the decoding accuracy was consistent within the F-zone, irrespective of navigational intent or visual cues during foraging (inbound vs. outbound). The spatial resolution encoded by the ensemble of neurons was precise enough to provide navigational utility (13.86 cm in the F-zone) considering the body length of an adult rat (∼20 cm excluding the tail). Therefore, we argue that the functional properties of mPFC neurons in the experiment are consistent with the spatial encoding hypothesis.

However, our observations also indicate that the spatial encoding in the mPFC was susceptible to disruption. Our distance regressor showed a higher error rate when rats engaged with the Lobsterbot inside the E-zone. A previous study also found that the presentation of a disfavored reward results in erroneous excursion of the location decoder near the feeder location (Mashhoori et al., 2018). Importantly, such disruption was not limited to the E-zone. Non-navigational behaviors such as grooming and sniffing greatly reduced the accuracy of location prediction suggesting that spatial encoding in the mPFC is tightly bound with task-demands such as goal-seeking and threat avoidance.

In the F-zone, where the regressor exhibited high spatial decoding accuracy, strategic navigation was primarily observed. Hok et al. (2005) argued that the mPFC neurons display spatial-selective firing patterns when the rat is engaged in navigation tasks rather than aimlessly wandering. Similarly, Ma et al. (2023) found that complex tasks involving working memory and reward improved spatial decoding accuracy in the mPFC. Consistent with these results, in our study, the removal of the Lobsterbot —which simplified the task to mere reward collection— resulted in decreased decoding accuracy. Given these observations, along with the mPFC’s lack of consistency in spatial encoding, it is plausible that the mPFC operates in multiple functional modes, and the spatial encoding mode is preempted when the complexity of the task requires deliberate spatial navigation.

Consistent with previous findings that mPFC neurons are sensitive to changes in valence such as reward (Horst & Laubach, 2013), threat (Burgos-Robles et al., 2009), or both (Kim et al., 2018), the recorded neurons modulated their firing rate around the entry and withdrawal to the E-zone, where sucrose reward and robotic threat coexisted. Consistent with the results indicating the mPFC is highly involved with executive functions (Dalley et al., 2004; G. W. Diehl & Redish, 2023), an ensemble classifier trained with the recorded unit activity predicted the success of AW at a fairly high accuracy (0.73). The prediction was significantly above the chance level (0.5), even several seconds before the action.

Interestingly, we found no difference between the PL and IL in the decoding accuracy of distance or avoidance decision. This somewhat surprising considering distinct roles of these regions in the long line of fear conditioning and extinction studies, where the PL has been linked to fear expression and the IL to fear extinction learning (Burgos-Robles et al., 2009; Dejean et al., 2016; Kim et al., 2013; Quirk et al., 2006; Sierra-Mercado et al., 2011; Vidal-Gonzalez et al., 2006). On the other hand, more Type 2 neurons were found in the PL and more Type 1 neurons were found in the IL. To recap, typical Type 1 neurons increased the activity briefly after the head entry and then remained inhibited, while Type 2 neurons showed a burst of activity during head entry and sustained increased activity. One study employing context-dependent fear discrimination task (Kim et al., 2013) also identified two distinct types of PL units: short-latency CS-responsive units, which increased firing during the initial 150 ms of tone presentation, and persistently firing units, which maintained firing for up to 30 seconds. Given the temporal dynamics of Type 2 neurons, it is possible that our unsupervised clustering method may have merged the two types of neurons found in Kim et al.’s study.

While we did not observe decreased IL activity during dynamic foraging, prior studies have shown that IL excitability decreases after fear conditioning (Santini et al., 2008), and increased IL activity is necessary for fear extinction learning. In our paradigm, extinction learning was unlikely, as the threat persisted throughout the experiment. Future studies with direct manipulation of these subpopulations, particularly examining head withdrawal timing after such interventions, could provide insight into how these subpopulations guide behavior.

We also found that the activity of Type 2 neurons was highly correlated with the success of AW. What is the role of the Type 2 sub-population in the execution of active defensive behavior? One straightforward hypothesis is that the increased activity simply represents internal motivation of fear or anxiety. This premise is supported by multiple experiments where the PL is vital for the expression of freezing responses (Burgos-Robles et al., 2009; Kim et al., 2018). Consequently, Type 2 units may represent the perceived threat level, which would strongly correlate with the rat’s voluntary head-withdrawal. Alternatively, the Type 2 sub-population may contribute to the computation or decision-making processes underlying defensive behaviors. This is based on a study by Halladay and Blair (2015), where they identified two sub-populations selectively responding to specific defensive behaviors, namely, conditioned movement inhibition and excitation. Because Type 2 neurons increased their activity in AW trials, it is possible that they overlap with the excitation-related neuronal subgroups in Halladay and Blair’s study. Despite this interesting conjecture, any analysis based on recording data is only correlational, mandating further studies with direct manipulation of the subpopulation to confirm its functional specificity.

Encoded representations in the mPFC range from prediction of upcoming decisions (Passecker et al., 2019) to rule switching (Laskowski et al., 2016). This multiplicity of cognitive functions in the mPFC has also been observed in primates, giving rise to the adaptive coding model (Duncan, 2001). This model posits that the mPFC neurons are fine-tuned to encode task-relevant concepts within specific contexts to meet the requirements of the current task. In our study, more than 80% of units—which coded distance information in the F-zone—modulated with the defensive decision making when entering or leaving the E-zone, thus supporting the adaptive coding model and functional switching between different contexts. The PCA results further revealed the heterogeneity of the populational activity in the E-zone as its manifold, or subspace of possible neural activity (for review, see Ebitz & Hayden, 2021), was distinctive from those of the F- and N-zones.

It is worth emphasizing that upon exiting the E-zone the neurons resumed encoding distance in the same manner as they did in previous trials. This is substantiated by the fact that our regressor was trained by the recording data from the whole session and therefore, it would have been impossible to accurately decode the distance if the neurons changed their encoding every time the rats re-entered the F-zone. Interestingly, in our foraging task, the rats quickly shifted their behavior mode between contexts—navigational in the F-zone vs. confrontational in the E-zone—within the same session. The temporal stability of spatial decoding in the F-zone throughout the session, despite the intermittent engagement in confrontational behavior in the E-zone, indicates that mPFC neurons selectively switch functions, rather than being randomly activated by the context change.

In summary, our findings reveal a dual, context-sensitive encoding in the rat’s mPFC, accommodating both spatial navigation and defensive decision-making during dynamic foraging. We demonstrated that the mPFC co-represented spatial and active foraging-related information depending on the context (Figure 8A) Figure. This functional switch not only highlights the mPFC’s multifaceted encoding capabilities but also reconciles conflicting models regarding its role in spatial orientation.

Hypothetical control models by which mPFC neurons assume different functional states.

(A) In the F-zone, navigational behaviors enhance the mPFC’s encoding of spatial information compared to other zones. In the N-zone, spatial coding diminishes when the rat engages in non-navigational behaviors. However, in the E-zone, these neurons shift their encoding strategy and become involved in coding for active foraging. We did not find a subset of neurons dedicated exclusively to either spatial coding or active foraging throughout the session. Instead, neurons changed their encoding scheme in a population-wide manner.

(B) Two hypotheses about how the switch is manifested. In this example, most mPFC neurons encode spatial information (blue circles). Information encoded in the mPFC can be regulated by internal/external arbitration signal (top-bottom blue arrow from green circles), or influenced by direct sensory inputs and navigation-related signals (left-right blue arrow) that prompt mPFC neurons to encode spatial information.

In primates and rodents, the mPFC is often likened to a “department store” due to its involvement in numerous cognitive functions (e.g., reward, decision-making, prediction error). One pioneering theory suggests that the mPFC’s involvement in diverse functions stems from its ability to learn associations between context, locations, events, and adaptive responses (Euston et al., 2012). However, in our multi-task design, decoder accuracy varied with context and task demands, prompting the question: how does the mPFC decide which information or associations to encode?

Two main hypotheses could explain this switch. A bottom-up hypothesis suggests sensory inputs or upstream signals dictate encoding priorities, while a top-down hypothesis proposes that an internal or external “arbitrator” selects the encoding mode and coordinates the relevant information (Figure 8B). Although the current study is only a first step toward finding the regulatory mechanism behind this switch, our control experiment, where rats reverted to a simple shuttling task, provide evidence that might favor the top-down hypothesis. The absence of the Lobsterbot degraded spatial encoding rather than enhancing it, indicating that simply reducing the task demand is not sufficient to activate one particular type of encoding mode over another. The arbitrator hypothesis asserts that the mPFC neurons are called on to encode heterogenous information when the task demand is high and requires behavioral coordination beyond automatic, stimulus-driven execution. Future studies incorporating multiple simultaneous tasks and carefully controlling contextual variables could help determine whether these functional shifts are governed by top-down processes involving specific neural arbitrators or by bottom-up signals.

Materials and Methods

Subjects

Eight male Sprague Dawley rats (250-300 g, 5-6 weeks of age, Orient Bio.) were used in this study. The rats were housed in pairs on a 12-h light/dark cycle with food and water available ad libitum for a one-week acclimation period. In the following week, rats were handled for 10 min each day until they displayed reduced fear responses towards the experimenter.

Apparatus and procedures

The apparatus followed the same dimensions as in Kimm & Choi (2018), measuring 800 mm × 480 mm and a height of 500 mm excluding the robot compartment. It was covered in gray EMF-shield fabric to reduce high-frequency noise during experiments. Battery-powered LED strip lighting illuminated the apparatus and −60 dB white noise masked background sounds. The system employed four Arduino Uno boards and a Raspberry Pi board for experiment control and behavioral tracking. All sensor and control signals were relayed to an RV5 multiprocessor (Tucker-Davis Technologies, FL, USA) for further analysis. Bird’s-eye view videos were captured during sessions, and the rat’s location was tracked using the butter package (https://github.com/knowblesse/butter). All control schematics and codes are publicly available (https://github.com/knowblesse/Lobster/Apparatus).

Before the training, food and water were gradually restricted until the rats’ weight reached 80% of their initial value. Over three days, they acclimated to the apparatus for 20 min. Training began with the habituation phase, during which rats had unrestricted sucrose port access for 20 min. The sucrose port has an IR sensor which was activated by a single lick. The rat usually stays in front of the lick port and continuously lick up to a rate of 6.3 times per second to obtain sucrose. Any sucrose droplets dropped in the bottom sink were immediately removed by negative pressure so that the rat’s behavior was focused on the licking. Upon reaching 1,500 licks per trial, rats proceeded to the three-day shuttling phase. The shuttling session was 30 min long, and the gate was closed 6 s after the first lick, forcing rats into the N-zone to reopen it. The subsequent Lobsterbot phase introduced the Lobsterbot, which executed two grabbing motions after 3 s (30%) or 6 s (70%) from the initial lick. After three days of Lobsterbot sessions, the rats underwent microdrive implant surgery, and recording data were collected from subsequent sessions, either Lobsterbot or shuttling sessions, depending on the experiment. For all post-surgery sessions, those with fewer than 20 approaches in 30 minutes were excluded from further analysis.

Surgery

Before surgery, rats were anesthetized using isoflurane (3%-5%) and oxygen (500 mL/min). Enrofloxacin (5 mg/kg) and meloxicam (1 mg/kg) were administered to prevent infection and mitigate postsurgical pain. A craniotomy was performed at AP +3.2 mm and ML 0.5 mm from the bregma using a carbide burr (Komet, Kempten, Germany) mounted on a high-speed drill (Strong 207A; Saeshin, Daegu, Korea). The tip of each microdrive was lowered to the following coordinates: AP +3.2 mm, ML right 0.3 ∼ 0.6 mm, DV −4 mm from the bregma. The precise ML coordinate was determined after the craniotomy to avoid damage to the superior sagittal sinus. Six to seven screws were securely affixed to the skull, and two additional screws above the cerebellum were used as the ground and reference points. Finally, dental cement was applied around the microdrive to secure the implant. Rats were allowed at least one week of recovery before the recordings. All procedures were approved by the Korea University Institutional Animal Care & Use Committee.

Data collection

Two types of microdrives were used in the recordings: a custom-made 4-tetrode microdrive and a pre-assembled 16-channel silicon probe (A2×2-tet-10 mm-150-150-121-HZ16_21 mm, NeuroNexus, MI, USA) mounted on a 3D printed microdrive. The microdrive featured an M1.6-18 mm screw with a 0.35 mm pitch. The recording site was retroactively determined by calculating the number of rotations from the final recording site, identified by a marking lesion.

During each recording session, the rat’s microdrive was connected to a ZIF-Clip digital headstage (Tucker-Davis Technologies, FL, USA). The neural signals were collected at a 25 kHz sampling rate, amplified, digitized from the head stage, and transferred to the RZ5 (Tucker-Davis Technologies, FL, USA) for data acquisition. The acquired data stream was further bandpass filtered between 300 and 5000 kHz, and all spikes above 3 to 4 sigma compared to the stream were collected for further clustering analysis. The Offline Sorter (Plexon, TX, USA) was used for the manual clustering process. All units with firing rates below 0.5 Hz or exactly 60 Hz (power noise) were considered artifacts and removed from further analysis. Sessions with fewer than 10 recorded units were also excluded.

Spatial decoding data preparation

For the neural data used in the spatial encoding analysis of the mPFC, the following method was employed. First, a Gaussian kernel with a standard deviation of 100 ms and a length of 1000 ms was prepared. Each unit’s spike data were down-sampled and serialized to a 1000 Hz vector, and the kernel was applied to smooth the signal. The mean and standard deviation values of this vector were computed and used to normalize the vector. After normalization, the vector was binned into 50 ms windows, resulting in 20 data points per second. Finally, each datapoint was paired with the corresponding location data.

The spatial firing map was constructed using the same normalized activity data. Neural activity was averaged for each pixel and then smoothed using a Gaussian filter with a sigma of 15 px and a filter size of 1001 px.

ANN regressor

A neural network regressor was built using the PyTorch package. It comprised seven layers in total, including the output layer; four of these were fully connected (FC) layers and two were drop-out layers. The initial weights for the FC layers were normalized with a mean of zero and a standard deviation of 0.2. Mean squared error loss was used in conjunction with the stochastic gradient descent method. For each regressor, a control regressor was built, which used the same data-label set, but the relationships were randomly shuffled. The evaluation of the regression results used a stratified 5-fold cross-validation while ensuring that each subset maintained the same proportion of data from the three zones. The five regression results per dataset were combined, and these results were compared with the true dataset using MAE. For the training/testing dataset, preprocessed 50 ms-long neural data were used. The binned neural data from all units (n) were concatenated to form a 1 × n size vector. These vectors were paired with the corresponding distance during the 50 ms time window. In the dataset adjusted for uneven location visits, we divided distance values into five equally sized bins. Then, a sub-dataset was created that contains an equal number of data points for each of these bins. Two different GPU servers were employed for training: a Tesla V100 from the National IT Industry Promotion Agency (NIPA, Seoul, Korea) and an in-house A4000.

Decoding error heatmap construction

To construct the decoding error heatmap, the following method was used. First, error data from each pixel were collected and mean values were computed. To accommodate pixels not visited by a rat during an entire session, an interpolation method was employed. A custom Matlab script, built using the scatteredInterpolant function, was used to interpolate errors for unknown locations. The 2D interpolated decoding error was subsequently smoothed using a Gaussian filter with a sigma value of 15 px and a filter size of 1001 px. These smoothed error values were then converted to cm, and a contour plot was overlaid using 25 levels.

PCA

A custom Python script for PCA was created using the sci-kit-learn package. Neural data vectors generated for the spatial decoding analysis were used for this analysis. The behavior and head tracking data were matched with neural data vectors, and the corresponding zone was labeled to the vectors. All neural data vectors from a single session were projected onto a new space using the PCA. For each zone, the centroid of the projected neural data vectors was calculated using Euclidean distance. The within-zone distance and the between-zone distance were calculated using the equations below.

where CZ denotes the centroid of the zone z and xZ represents a neural data vector from the zone z.

Representative figures were drawn in a manner to highlight the cluster differences. Among all projected neural vectors designated to a zone, the 2000 vectors with the shortest distance from the centroid were selected. From these vectors, 200 were randomly chosen and the first two dimensions of the projected value were plotted.

Behavior-responsive unit analysis

Behavior-responsive units were identified using the following method. For each trial, spikes within a 4-s window (2 s before and after the onset of three events: gate opening upon returning to the N-zone, head-entry, and head-withdrawal) were counted and registered in 50 ms-long time bins, resulting in 80 data points per trial. These neural activities were then averaged across all trials, producing a 1 × 80 vector for each event. The binned neural activity ± 2 s from the entrance to the N-zone was chosen as a baseline, considering that the rat was consistently oriented towards areas other than the robot. The mean and standard deviation of the baseline were used to normalize activity during the head-entry and the head-withdrawal. The overall method created two event vectors, head-entry-vector and head-withdrawal-vector, per unit. A unit was considered behavior-responsive to an event if any value within the 80 data points of the event vector was greater or less than 3 sigma.

Hierarchical cell clustering analysis

To cluster cells that exhibited similar responses during the head-entry and/or head-withdrawal, an unsupervised hierarchical cell clustering method was employed. For every pair of units, the similarity score using Pearson’s r was calculated, and the pair with the highest similarity score was grouped to form a new pseudo-unit that had the average response of the two units. This procedure was repeated until all units were grouped into eight pseudo-units. Among each pseudo-unit, if the number of units belonging to the group was equal to or greater than 50, the pseudo-unit was considered a valid cluster. Otherwise, the pseudo-unit was classified as an artifact.

To compensate for the limited number of neurons recorded per session, the hyperparameter set was chosen to generalize their activity and categorize them into major types, allowing us to focus on neurons that appeared across multiple recording sessions. Although changes in the hyperparameter sets resulted in different numbers of clusters, the major activity types remained consistent (Supplementary Figure S8). However, there is a chance that this method may not differentiate smaller subsets of neurons, particularly those with fewer than 50 recorded neurons.

Run-and-stop event analysis

To examine whether a change in movement velocity was related to neural activity, a “run-and-stop” event was defined by the following criteria: 1) to qualify as “run” at a given moment, the rat must be moving at least 7.5 cm/s for at least 2 s or longer; 2) to qualify as the “stop” at a given moment, the rat must be moving at most 5 cm/s for at least 1 s or longer; and 3) the “run-and-stop” event was the time point between “run” and “stop” outside the E-zone. The temporal relationship between the run-and-stop event and neural activity was visualized using the same method as in the behavior-responsive unit analysis.

Naïve Bayes AW/EW classifier and data preparation

For the neural data used in the head-withdrawal-type (AW/EW) classification analysis, a data preparation method similar to that of the spatial decoding analysis was employed. First, the same Gaussian kernel was generated and applied to the down-sampled and serialized neural data. Mean and standard deviation values of this vector were computed for later normalization. Unlike the 50 ms-length neural data used in spatial regression, 2-s-length neural data was used. For each trial, all spikes within a 2-s window around the onset of the head-withdrawal were collected and serialized into a 1000 Hz vector. The Gaussian kernel was applied and normalized using precalculated values. This vector was then binned into 50 ms intervals, resulting in 40 data points per event. Binned vectors from all recorded units (n) were concatenated to form a 1 × (40n) size vector per trial. These vectors were then paired with their corresponding head-withdrawal-type, AW or EW. In the analysis using varying time windows, the previously employed data preparation method was the same, but the time range of the neural data varied relative to the head-withdrawal.

A naïve Bayesian classifier based on the multivariate Bernoulli model was built using the sci-kit-learn package. A custom Python script was written to train and evaluate the classifiers. The classifier employed a uniform prior and smoothing parameter, alpha, which was set to 1 in all classifiers used in the study. Similar to the analysis using the ANN regressor, a control classifier was built that was trained with shuffled data. The evaluation of the classifier results employed the stratified 5-fold cross-validation method with each subset maintaining the same proportion of head-entry/head-withdrawal. Five classification results per dataset were combined, and these results were compared with the true dataset using balanced accuracy, calculated using the equation below.

To interpret how the classifier anticipated upcoming head-withdrawal behavior, class discrimination index was employed, focusing on a dataset with a time window of −2 to 0 s. To calculate the index, the class likelihood for every dimension was extracted during each test dataset evaluation, and a natural logarithm was applied to balance the small likelihood values. With each unit’s neural data divided into 40 bins, the log likelihood values were added to create a joint probability. This process yielded five joint likelihood values per unit, and their average was computed. Following Mladenić & Grobelnik’s (1999) method, the log odds ratio between P(unit|AW) and P(unit|EW) was calculated for all units using the equation below.

where P(unit|AW) denotes joint likelihood (probability of observing unit activity given the AW) and P(unit|EW) denotes the same for the EW. These values were then compared according to the unit’s assigned group: Type 1, Type 2, or others. Taking the natural log of the odds ratio brings about a symmetricity in the results. A positive value represented a unit’s relatively increased activity during the AW compared to the EW, and a negative value represented the opposite.

Statistical analysis

Statistical analyses were performed using Prism (GraphPad, MA, USA). To confirm normality of the data, Shapiro-Wilk test was applied to all data we obtained from the experiment. For all statistical tests, p-value of <.05 were used to assess the significance of differences and all t-tests were two tailed. The correlation analysis comparing feature importance scores from two decoders excluded the 1% of outliers before computing Pearson’s r. Both ANOVA and mixed-effect model analysis applied the Geisser-Greenhouse correction to account for violations of the sphericity assumption. The corrected degree of freedom is presented alongside F values. For multiple comparison tests, Sidak’s correction was used to control the family-wise error rate. For multiple linear regression analysis of the distance decoding error, the following model was incorporated.

The model fitting method incorporated the least squares method assuming a Gaussian distribution of residuals. Boolean variables, presence of Lobsterbot, and recording site used dummy codes (0: Lobster, 1: Control; 0: PL, 1: IL) to represent the state of the experiment. All error bars in graphs represent one standard error of the mean. Statistically significant differences are denoted as * for p <.05, ** for p <.01, and *** for p <.001.

Supplementary Figures

Head withdrawal time distribution across all subjects, categorized by trial type

Despite the use of two distinct attack times (3s and 6s), there was no noticeable increase in head withdrawals around the 3-second mark, indicating that the rats did not rely on a 3-second cou ntdown. Instead, they exhibited a relatively stable distribution leading up to the 6-second attac k.

Foraging-related behavioral indices fluctuate upon the initial encoun ter with the Lobsterbot but stabilize after 3 sessions.

(A) Number of approaches. The number of approaches, measured in total trials, decreased afte r the initial encounter with the Lobster (Lob1), but later increased after 3 Lob sessions. (B) Nu mber of licking behaviors. The number of licking behaviors significantly decreased during the first encounter but returned after 3 sessions. (C) Number of licks per trial. The number of licki ng behaviors per trial was decreased after the encounter. (D) Lick latency. The lick latency incr eased after encounter but returned to pre-encounter level after 3 sessions. The black dotted line indicates the timepoint of surgery. Error bars represent SEM.

Comparison between distance regressor algorithms

A comparison of distance regressors using the same dataset showed that the artificial neural network (ANN) and the random forest regressor outperformed the others, but ANN was chosen for its strong generalization to noisy neural data and robustness to hyperparameters. Error bars represent SEM.

Distances between each centroid pairs from all recording sessions. Figure 4B result was reanalyzed using partial dataset which excluded “critical event times” defi ned as ±1 second from head-entry and head-withdrawal. Even removing behaviorally significant data, E-zone’s populational activity was distinctly positioned compare to other zones. Error bars repres ent SEM.

A run-and-stop event (sudden velocity drop outside the E-zone) does not evoke neural modulation.

Normalized activity of HE1 and HE2 units during run-and-stop events (colored; HE1-r&s and HE2-r&s) show no modulation of neural activity compared to highly modulated activity around the head-entry (black and gray; HE1-HE and HE2-HE). Gray shadings indicate SEM.

Most units are classified into either the HE1-HW1 or HE2-HW2 gro ups

(A) Confusion matrix comparing the Head Entry and Head Withdrawal groups. A large propo rtion of units fall into either the HE1-HW1 category (n=299) or the HE2-HW2 category (n=94).

(B) Normalized neural activity of Type 1 (HE1-HW1) and Type 2 (HE2-HW2) neurons during the head-entry and the head-withdrawal. Gray shadings represent SEM.

Type1 and Type 2 neurons’ PETHs around head withdrawal separated by AW and EW

(A) (Top) Type1 neurons’ PETH around head withdrawal separated by AW and EW. Each PETH is sorted by neuron’s peak timepoint. (Bottom) Average of PETH. (B) Same as (A) with Type 2 neurons.

Hierarchical clustering results with different hyperparameter sets Hierarchical clustering uses two hyperparameters: cutoff limit and the number of initial clusters. These variables were varied to assess the clustering results’ dependence on them.

While changes in these variables affected the number of groups, the response characteristics of the top two groups (which were used to further classify Type 1 and Type2) remained consistent. Gray shadings represent SEM.

AW/EW prediction accuracy as a function of dataset time window.

Data Availability

All raw data including behavior, neural, and video are available upon request.

Acknowledgements

This work was supported by the National Research Foundation of Korea (NRF-2021M3E5D2A01023887) funded by the Korean Government, MSIT. GPU resource was supported by High-performance computing support project by NIPA, Seoul, Korea.

Additional information

Code availability

The Matlab and Python scripts used for classification and behavioral analysis are openly accessible at the following repository: https://github.com/knowblesse/Lobster.

Author Contributions

J. H. J and J-S. C designed the study. J. H. J conducted the experiment and analyzed the data. J. H. J and J-S. C wrote the manuscript.

Additional files

Video 1

Video 2

Video 3