Population analyses reveal heterogenous encoding in the medial prefrontal cortex during naturalistic foraging

Foraging in the wild requires coordinated switching of critical functions, including goal-oriented navigation and context-appropriate action selection. Nevertheless, few studies have examined how different functions are represented in the brain during naturalistic foraging. To address this question, we recorded multiple single-unit activities from the medial prefrontal cortex (mPFC) of rats seeking a sucrose reward in the presence of a robotic predator (Lobsterbot) that posed periodic threats. Simultaneously recorded ensemble activities from 10-24 neurons were analyzed in reference to various behavioral indices as the animal moved freely between the nest (N) and the goal (E) across the foraging (F) area. In the E-zone, the rat initially received and gradually learned to avoid unpredictable attacks by the Lobsterbot. An artificial neural network, trained with simultaneously recorded neural activity, estimated the rat’s current distance from the Lobsterbot. The accuracy of distance estimation was the highest in the middle F-zone in which the dominant behavior was active navigation. The spatial encoding persisted in the N-zone when non-navigational behaviors such as grooming, rearing, and sniffing were excluded. In contrast, the accuracy decreased as the animal approached the E-zone, when the activity of the same neuronal ensembles was more correlated with dynamic decision-making between food procurement and Lobsterbot evasion. A population-wide analysis confirmed a highly heterogeneous encoding by the region. To further assess the decision-related activity in the E-zone, a naïve Bayesian classifier was trained to predict the success and failure of avoidance behavior. The classifier predicted the avoidance outcome as much as 6 s before the head withdrawal. In addition, two sub-populations of recorded units with distinct temporal dynamics contributed differently to the prediction. These findings suggest that the mPFC neurons may adopt at least two modes of heterogenous encoding that reflect the processing of relevant spatial context and the imminent situational challenge.


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 taskswitching (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., 2007Mobbs et al., , 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 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 populationlevel 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.

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

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.
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 machinelearning 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 50ms 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 ttest, 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.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 PL and IL might encode the spatial distance from the goal.(D) Prediction accuracy in the Fzone 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.

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 Ezone, 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).
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 nonnavigational 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 nonnavigational 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.
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 Ezone 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.
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 subpopulation (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).
Based on the previous clustering analysis that identified functionally cohesive subpopulations 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 subpopulations 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 taskrelevant 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 Ezone-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.In the F-zone, goal-directed navigation prompts a greater number of mPFC neurons to encode spatial information compared to other zones.In the E-zone, these neurons switch the encoding scheme for reward-or threat-related information.

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.Batterypowered 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 highspeed 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 subdataset 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   denotes the centroid of the zone  and   represents a neural data vector from the zone .
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, HEvector 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.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.

Figure 1 .
Figure 1.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 a AW trial, the rat typically retracts its head ahead of time and watches the Lobsterbot attack.In a 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 L sessions.

Figure 2 .
Figure 2. Ensemble activity from PL and IL predicts distance from the goal.(A) Schematic diagram of electrode implantation and estimated recording site.Top: A movable 4tetrode 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) 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 Lobsterbot.The example data depicted in the figure is a sample recording from 20 units when the rat is at a particular distance away from Lobsterbot, indicated by the white bold line.(C) Accuracy of the location prediction.Mean Absolute Error (MAE) was calculated for the two types of decoders: one with the original dataset (Original) and

Figure 3 .
Figure 3. Spatial encoding is disrupted by non-navigational behaviors (NNBs).(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 Nand 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 NNBs 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 location-decoding errors (N-zone) during navigational vs. NNBs.The error was significantly larger when the rat was engaged in NNBs within the N-zone.(E) Comparison of location-decoding errors of regressors trained with and without NNBs.The overall decoding error was significantly smaller when the ANN was trained without the data during NNBs.

Figure 4 .
Figure 4. 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.

Figure 5 .
Figure 5. Peri-Event Time Histogram (PETH) of all units -responsive to the head entry (HE), aligned by the peak firing rate and the averaged activity.(A) Top: The PETH for each unit is based on the Z-score of which magnitude is color-coded (temperature bar on the right).Bottom: The red vertical lines mark the timing of the HE.The peak latency of each unit varies from as early as 2 s before and to 1~2 s after the HE.(B) Functional segregation of the 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) PETH of all units responsive to the head withdrawal (HW) and the averaged activity.(D) Functional segregation of the recorded units.Top and middle: Three sub-populations of units based on hierarchical cell clustering analysis.Bottom: The averaged activity for each sub-population.

Figure 6 .
Figure 6.Neural ensemble activity predicts failure and success of avoidance response.(A) Schematics of the event decoding analysis.The Naïve Bayesian decoders (NBCs) 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 NBC.The decoding accuracy of the classifier was significantly higher than that from the shuffled data.(C) Temporal characteristics of prediction accuracy of the NBC.Prediction accuracy was significantly higher at the time points as early as 5-7 s before the HW.(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.

Figure 7 .
Figure 7. mPFC neurons switch between different functional modes depending on the location.
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 analysisStatistical 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.MSE = β 0 + β 1 (number of units used in the decoder) + β 2 (presence of Lobsterbot) + β 3 (recording site; PL, IL)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.