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 plays 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 deliberately strategizing defensive behaviors across different situations (see Mobbs et al., 2020 for review), navigation during naturalistic foraging where risk and reward coexist may serve as a suitable context for exploring the flexible information processing capabilities of the mPFC. In fact, the mPFC houses spatially selective neurons that represent relative maze positions (Kaefer et al., 2020; Yang & Mailman, 2018), goal locations (Hok et al., 2005; however see Böhm & Lee, 2020 for contrary results), and precise Cartesian coordinates during food procurement (Mashhoori et al., 2018).

The mPFC plays a role in processing emotional stimuli, especially those signaling threats in the environment. Not only is the mPFC intricately intertwined with the major limbic areas including the amygdala (Gabbott et al., 2005; Hoover & Vertes, 2007; Vertes, 2004), but it is also causally linked to defensive behaviors. For example, the infralimbic cortex (IL) of the mPFC is associated with the extinction of conditioned fear memory (Sierra-Mercado et al., 2006, 2011). 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). These findings suggest that the critical information for risky decision-making might be represented and processed in the mPFC. Furthermore, PL neurons show increased activity in response to survival-related stimuli such as food, artificial predators, or both (Kim et al., 2018).

Currently, there is limited understanding of how information is represented in the mPFC during dynamic foraging. It is unclear whether individual cells are dedicated to spatial representation and emotional stimuli processing or if they adapt their function to the current goal. In this study, rats engaged in a naturalistic foraging task and produced defensive responses to a robotic predator. This approach allowed continuous monitoring of the neural activity as the rat switched between navigation and strategic reward procurement. By employing population-level and single-level analyses, we explored the mPFC’s role across multiple behaviors and observed a functional mode switch characterized by the co-representation of different functions in distinct spatiotemporal regions.

Results

Mix of avoidance and escape behaviors during naturalistic foraging under threat

The foraging arena, modified from (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 walled passage (Figure 1A): the foraging (F) and nest (N) zones. At one end of the F-zone, a robot compartment was located, separated by a motorized gate. The part of the F-zone where a rat could push its body into the robot compartment was separately named the encounter zone (E). The entrance to the E-zone was monitored by the beam sensor mounted near the gate. Inside the robot compartment, a sucrose port was situated, activated by licking, and guarded by a robotic claw (called Lobsterbot owing to its striking attack; Figure 1B).

Figure 1.

The entire experiment was composed of three distinct phases (Figure 1C): habituation (H), shuttling (S), and foraging under Lobsterbot threat (L). The L phase was further divided into early (L1) and late (L2) phases separated by the time of surgery. During the H phase, a food-deprived rat was habituated to the sucrose port that released a drop of sucrose solution activated by licking. Once the rat maintained a stable licking frequency, typically after 3 days, the next phase (S) started. During the S phase, the rat was trained to shuttle back and forth to procure sucrose. The robot compartment’s gate was closed 6 s from the first lick and would only reopen once the rat returned to the N-zone. The S phase lasted 3 days, followed by the introduction of the L phase. During the L phase, the robot executed a fast-striking attack, 3 s (30%) or 6 s (70%) from the first lick, varying randomly each trial to mimic natural uncertainty and prevent simple time-counting strategies. As in the S 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 L phase, the foraging rat was given a chance to 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 the 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).

All five rats exhibited distinctive phasic behavioral changes. Notably, there was a significant drop in approach attempts (i.e., number of trials) during the first session of the L1 phase (Lob1), followed by a stable rate of approach attempts in subsequent sessions (Supplementary Figure S1A). Other behaviors also followed similar phasic changes, including the number of licks (Supplementary Figure S1B), licks per trial (Supplementary Figure S1C), and lick latency after gate opening (Supplementary Figure S1D). To assess these behavioral changes, we compared the last session of the S phase (Sht3) and all sessions of the L1 phase (Lob1, Lob2, Lob3). The analysis showed that the number of approach attempts (repeated measure one-way analysis of variance (ANOVA); F(1.29, 5.16) = 7.18, p =.038), the number of lickings (F(1.64, 6.56) = 12.87, p =.006), and the licks per trial (F(2.21, 8.84) = 23.13, p <.001) was significantly changed during the L1 phase, while the lick latency (mixed-effects analysis; F(1.11, 3.71) = 4.06, p = .119) only showed numerical change. More specifically, approach attempts and licking were significantly decreased at the first session of the L1 phase (Sht3 vs. Lob1; t(4) = 4.53, p = .031; t(4) = 5.81, p =.013) and subsequently increased at the end of the L1 phase (Lob1 vs. Lob3; t(4) = 8.13, p = .004; t(4) = 4.84, p =.025). In summary, after the three L sessions, the rats’ sucrose consumption reverted to S phase levels (Sht3 vs. Lob3; t(4) = 1.16, p = .670) by increasing their number of approaches, compensating for a numerically decreased number of licks per trial (Sht3 vs. Lob3; number of approaches: t(4) = 6.13, p = .011; number of licks per trial: t(4) = 3.46, p = .075).

In the early sessions of the L1 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 end of the L1 phase. Compare to the S phase, where the rats remained in the E-zone until the gate closed in 88.96% trials, they began to display a mix of AW and EW by the end of the L1 phase (AW ratio; 65.11%). Due to the transient nature of the decreased approach attempts and the various defensive behaviors after the first encounter, surgery was performed after three L sessions when the rats displayed a range of 30% to 60% AW (Figure 1F).

Population activity in the PL and IL predicts the distance from the goal while navigating 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 that was performed after the last 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.

Figure 2.

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 the mPFC units had spatial correlates, we performed a regression analysis for the entire arena (N-, F-, and E-zones). The inputs to the regressor were the simultaneously recorded neural activity and the predicted output was the distance from the goal location computed from the rat’s head position. Among machine-learning algorithms, we implemented a four-layer artificial neural network (ANN; see Materials and Methods for a detailed structure), as the ANN achieves significantly fewer decoding errors compared to other algorithms (Mashhoori et al., 2018). More specifically, to decode spatial information from neural data, recorded neural activity was segmented into 50-ms bins and used as input; the distance from the center of the robot was used as the output (Figure 2B).

The analysis of trained ANNs revealed that the distance from the rat to the robot 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), indicating that the population activity in the PL and the IL contain spatial information (Figure 2C, Video 2). 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. Owing to 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 spatial encoding 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 2D). In summary, these data suggest that precise spatial encoding could be extracted using the neural ensemble in the mPFC.

The decoding accuracy remained unchanged despite the differences in motivation and perceived visual cues due to the movement direction.

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 goal-oriented navigation behavior. 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 S phase, where the Lobsterbot was removed. 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 (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 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.

Excluding non-navigational behavior (stereotypic behavior) enhances spatial accuracy of the population decoding in the N-zone

Although the initial analysis showed that the population activity encoded the distance of the rat to the goal reasonably well, its 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).

Figure 3.

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).

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 spatial encoding in the N-zone might have resulted from the disruption induced by non-navigational behaviors.

Multi-dimensional 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-zone 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.

Figure 4.

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 one-way repeated measure ANOVA showed a significant difference between centroid pairs (F(1.61, 62.64) = 51.99, p <.001). Further post hoc analysis revealed a significant difference between N-F and N-E pairs (t(39) = 9.96, p <.001), N-F and F-E pairs (t(39) = 2.55, p =.044), and N-E and F-E pairs (t(39) = 9.13, p <.001), 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 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. In summary, populational activity inside the E-zone showed distinctive activity patterns unlike those of active navigation.

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

In the E-zone, navigational behavior was mostly replaced by approach attempts and avoidance. 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 (HE) and the head withdrawal (HW). 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 event-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 HE and HW.

Figure 5.

For the HE, 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 HE, 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 HW. 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 HW (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 (HE and HW), 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 for the HE (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. For the HW, 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 HW. The second (HW2: n = 147, 23.26%) maintained increased activity until the onset of HW 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 HW.

Well-trained rats typically run across the F-zone at high velocity and arrive at the sucrose port (HE) in the E-zone, stopping abruptly. Moreover, the HW 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 HE and/or HW. To exclude this possibility, a “run-and-stop” event was defined and correlated with the neural data. This 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 S2, there was no significant neural modulation detectable around the run-and-stop events. These findings suggest that the neural modulation occurring around the HE and HW is not attributable to changes in velocity of movement.

Interestingly, a significant categorical overlap between HE and HW neurons was found. In the HW1 group, 78.68% (n = 299) of units were classified as HE1. In the HW2 group, 63.95% (n = 94) of units were from HE2 (Supplementary Figure S3A). 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 S3B). Specifically, Type 1 neurons exhibited a rapid decline in activity following a short peak 200 ms before the HE, whereas Type 2 neurons displayed a sustained high peak activity (z ≍ 2) from the onset of the HE until the HW. 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 withdrawal-related activity to determine their functional contributions to avoidance decisions. If the recorded neurons in the mPFC encode avoidance decisions, AW and EW may be distinguishable from decoded population activity. A naïve Bayesian classifier (NBC, 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).

Figure 6.

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 NBC 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 NBCs using 2-s ensemble activity at varying time points before the execution of the HW 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 (multiple paired t-tests with Šídák-Bonerroni correction; -31 – -29: t(39) = .24, p > .999; -21 – -19: t(39) = .57, p > .999; -11 – -9: t(39) = 1.62, p = .701; -7 – -5: t(39) = 5.381, p < .001; -6 – -4: t(39) = 5.328, p < .001; -5 – -3: t(39) = 4.610, p < .001; -4 – -2: t(39) = 4.631, p < .001; -3 – -1: t(39) = 4.544, p = .008; -2 – 0: t(39) = 4.544, p < .001; -1 – +1: t(39) = 8.505, p < .001; Figure 6C).

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 NBC 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 HW (from 2 s before HW to the onset of HW) 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 HW behavior, and increased activity of unit subgroups was strongly associated with the execution of AW.

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 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 built an ANN regressor targeting the distance, a more generalized spatial index which also reflects an internal representation of value and risk.

The decoding analysis showed that brief navigation was sufficient to engage a majority of mPFC neurons recorded during foraging into spatial encoding. Furthermore, the spatial encoding within the F-zone showed a consistent representation, irrespective of navigational intent or visual cues during foraging (inbound vs. outbound). The spatial resolution encoded by the ensemble of neurons was precise (13.86 cm in the F-zone) considering the body length of the adult rat used in the experiment (20 cm or above excluding the tail). Therefore, we argue that mPFC neurons 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, appearing during non-navigational behaviors unrelated to reward or goal and significantly degrading the quality of spatial representation. These results suggest that spatial encoding in the mPFC was task-dependent.

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 some of which involve spatial encoding.

Consistent with mPFC neurons being responsive to 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. As the mPFC is highly correlated with executive functions (Dalley et al., 2004; G. W. Diehl & Redish, 2023), we built an ensemble classifier and found that mPFC activity predicted the success of AW, possibly by integrating sensory information represented inside the region.

We also found that the high activity levels of Type 2 neurons were 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 fear or anxiety motivation. 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 HW. 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.

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 in previous trials. This is substantiated by the fact that our regressor was trained on whole session data rather than trial-by-trial data, and it would have been impossible to accurately decode distance if the neurons changed their encoding method every time the rats re-entered the F-zone. Interestingly, in our experiment paradigm, the rats quickly shifted between contexts—navigational in the F-zone vs. confrontational in the E-zone—within the same session. The consistent spatial decoding in the F-zone throughout the session indicates that mPFC neurons selectively switch functions, rather than being randomly activated when the context changes.

In summary, our findings reveal a dual, context-sensitive functionality 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 foraging-related information—such as reward, threat, and execution of defensive behavior—depending on the context (Figure 7). This functional switch not only highlights the mPFC’s multifaceted encoding capabilities but also reconciles conflicting models regarding its role in spatial orientation. Future studies should be performed to determine whether these functional shifts are governed by top-down processes involving specific neural arbitrators, or by bottom-up processes driven by signal strength.

Figure 7.

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, daily confined from the robot compartment for 20 min. Training began with the H phase, during which rats had unrestricted sucrose port access for 20 min. Upon reaching 1,500 licks per trial, rats proceeded to the three-day S phase. The S phase 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 L phase introduced the Lobsterbot, which executed two grabbing motions after 3 s (30%) or 6 s (70%) from the initial lick. Lasting three days, rats attempting fewer than 30 approaches in 30 min were excluded from further experiments.

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 (KUIACUC-2020-0029).

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 (A2x2-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.

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 C^z denotes the centroid of the zone z and x^zrepresents 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, HE, and HW) 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 HE and HW. The overall method created two event vectors, HE-vector and HW-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 HE and/or HW, 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. Although this method had a chance of losing small population groups that exhibited important unique activity during the event, other groups were focused on as they consistently appeared across multiple recording sessions and subjects.

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 HW-type classifier and data preparation

For the neural data used in HW-type 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 HW 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 HW-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 HW.

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 HE/HW. 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 HW 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) software. All paired t-tests used two-tailed p-values to assess the significance of differences. A p-value of <.05 was considered statistically significant. 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.

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.

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.

Author Contributions

1. 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.

Competing Interest Statement

The authors declare no competing interests.

Supplementary Figures

Supplementary Figure 1

Supplementary Figure 2

Supplementary Figure 3

Video 1. https://public.knowblesse.com/eLife2023/Video1

Video 2. https://public.knowblesse.com/eLife2023/Video2