Authors: He Chen, Jun Kunimatsu, Tomomichi Oya, Yuri Imaizumi, Yukiko Hori, Masayuki Matsumoto, Yasuhiro Tsubo, Okihide Hikosaka, Takafumi Minamimoto, Yuji Naya, Hiroshi Yamada
Categories: Article, Biological sciences, Natural sciences, Neuroscience, Systems neuroscience
Source: iScience
Authors: He Chen, Jun Kunimatsu, Tomomichi Oya, Yuri Imaizumi, Yukiko Hori, Masayuki Matsumoto, Yasuhiro Tsubo, Okihide Hikosaka, Takafumi Minamimoto, Yuji Naya, Hiroshi Yamada
Neural dynamics are thought to reflect computations that relay and transform information in the brain. Previous studies have identified the neural population dynamics in many individual brain regions as a trajectory geometry, preserving a common computational motif. However, whether these populations share particular geometric patterns across brain-wide neural populations remains unclear. Here, by mapping neural dynamics widely across temporal/frontal/limbic regions in the cortical and subcortical structures of monkeys, we show that 10 neural populations, including 2,500 neurons, propagate visual item information in a stochastic manner. We found that visual inputs predominantly evoked rotational dynamics in the higher-order visual area, TE, and its downstream striatum tail, while curvy/straight dynamics appeared frequently downstream in the orbitofrontal/hippocampal network. These geometric changes were not deterministic but rather stochastic according to their respective emergence rates. Our meta-analysis results indicate that visual information propagates as a heterogeneous mixture of stochastic neural population signals in the brain.
Visual inputs activate a large number of neurons in the brain that construct numerous neural networks to process information in an environment.^1^^,^^2^^,^^3^ This brain-wide activity change reflects the information processing embedded in each individual neural circuit; however, limitations of spatial and temporal resolution in the measurements of circuitry activity disrupt our understanding of brain-wide visual information processing.^4^^,^^5^^,^^6^^,^^7^^,^^8^^,^^9^ Under this limitation, considerable attempts have been made toward understanding how the brain processes information using a variety of developing theoretical frameworks.^10^^,^^11^^,^^12^^,^^13^^,^^14^^,^^15^
One of the analytic frameworks developed within the last two decades is state-space analysis,^16^ which provides a mechanistic structure of information processed in the lower-dimensional space of a neural population.^17^^,^^18^^,^^19^ This analytical tool identifies dynamic neural population structures that reflect information processing for general biological features^20^^,^^21^ and allowed us to describe those features as a neural geometry with high temporal- resolution^13^^,^^14^^,^^15^ in the sub-second order. A large number of neural circuits may process information moment by moment,^6^ and they may form a population geometry, such as rotational,^18^ curvy^19^^,^^22^ or straight^23^ geometries as the typical and basic features of dynamics. A recent finding suggests that the combination of neural population geometries may be the key to processing information to transform sensory inputs into memory.^24^ Recent studies have extended these analytical frameworks^25^^,^^26^^,^^27^^,^^28^^,^^29^; however, because a limited number of studies compared among brain-wide neural populations, how the brain processes information in the form of geometry remains elusive (except^30^).
To examine how brain-wide neural dynamics are formed to process visual information, we accumulated the neural population data of monkeys from four laboratories that contained 10 neural populations, including 2,500 neurons across temporal/frontal/limbic networks (i.e., meta-analysis). We applied the state space analysis (i.e., dimensionality reduction) to the neural population data followed by a bootstrap resampling technique that detects and replicates neural modulation dynamics in a low-dimensional neural space as a parental population. Following the analysis of bootstrap replicates with the Lissajous curve function, our cross-study comparison revealed that a gradual shift in stochastic neural population signals occurred throughout the temporal-to-frontal brain regions.
We compared the trajectory geometries across many neural populations widely distributed in the brain from the output brain regions in the ventral visual pathway^31^^,^^32^^,^^33^^,^^34^ and its downstream brain regions that may access memories associated with a visual stimulus. These included ten brain regions that were accumulated from nine monkeys examined in the four laboratories (Table S1), from the higher-order visual area TE and their downstream brain regions in cortical, subcortical, and limbic structures, such as the temporal/orbitofrontal cortices, striatum, and hippocampus (HPC) (Figure 1A, No. 1 to 10). A total of 2,500 neurons were accumulated across the four behavioral tasks (Figure S1), in which visual items provided monkeys with position and/or reward information during the active (Exps. 1 and 3) and passive (Exps. 2 and 4) behavioral responses. All these behavioral tasks require monkeys to perceive visual cues to perform the Exp. 1. Item-location-retention (ILR) task; Exp. 2. Scene-based object-value task; Exp. 3. Delayed reward tasks; Exp. 4. Single cue task (see STAR Methods). Using state-space analysis, we characterized structures of neural population geometry that appeared in the lower dimensional neural space, which describes how neural modulation by the task parameters of interest processes information at the population level^23^^,^^29^ (i.e., targeted dimensionality reduction^35^).Figure 1Neural population geometries in the visual memory pathway(A) Anatomical depiction of neural populations obtained from the 10 brain regions in nine macaques during the four different behavioral tasks in Exps. 1 to 4.(B–E) Rotational (B), curvy (C), straight (D), and unclear dynamics (E) detected by visual inspection. Single trajectory geometry was obtained in each brain region using neural population data. Number of neurons were 116–590 (Table S1). In A-E, the 10 brain regions are numbered as 1. TE, 2. STRt, 3. PRC, 4. CDb, 5. HPC, 6. VS, 7. cOFC, 8. CDh&b, 9. PHC, and 10. mOFC. The 0.05 s time bin was used for the analysis. The time from visual stimulus appearance was shown in sec. Characters indicate stimulus conditions; Ib: best item; I2 to I7: second best to seventh best item; Iw: worst item; D: delay, M: magnitude; P: probability. See also Figures S1–S3 and Table S1.
We found all three types of geometric patterns, rotational, curvy, and straight geometries, in the top two dimensions (Figures 1B–1D, see Figure S2 for performance as the percentage of variance explained), including unclear structures based on visual inspection (Figure 1E). All 10 neural populations showed a significant structure at the principal component (PC)1–2 plane based on shuffle controls (Figure S3, p < 0.05 for all PC1 and PC2). PC3 did not show statistical significance in some neuronal populations (Figure S3, Exps. 2 and 4; compare the black and gray dots for each PC). These identified geometric structures appeared to be distributed from complex to simple, reflecting the circuit distance to the visual input (Figure 1). For example, TE (Figure 1B, No. 1, A11 plane in A) and its downstream region striatum tail (STRt) (Figure 1B, No. 2, A11 plane in A), showed rotational geometries (i.e., circle) during visual item recognition. In more detail, at the beginning of information processing after the visual item presentation (Figure 1B left, see green s at time = 0), the neural state was positioned at around the center of the PC1-2 plane, and then rotated at approximately 0.2 s and started to move the second quadrant with the counterclockwise rotation going back close to the initial point (see green e 0.6 s after visual item onset). An opposite rotation was observed for the worst visual item, with a smaller change (orange). These rotational structures were also observed at the STRt on a similar timescale (Figure 1B, right).
In the downstream brain regions, such as the perirhinal cortex (PRC) and caudate body (CDb), rotational or curvy dynamics were observed (Figure 1C, No. 3, PRC, A11 plane in A and No. 4, CDb, A23 plane in A), which were characterized by a half rotation ending at the opposite neural space and endpoints deviating from the initial point (i.e., half cycle), although it was unclear whether they showed rotational or curvy dynamics. In contrast, straight dynamics were observed in brain regions far away from the visual inputs in the HPC (Figure 1A, No. 5 in the A11 plane; Figure 1D) and the central part of the orbitofrontal cortex (cOFC) (Figure 1A, No. 7 in the A32 plane; Figure 1D). In addition, the ventral striatum (VS) showed straight dynamics (Figure 1A, No. 6 in the A23 plane; Figure 1D), although some structures could not be clearly determined (Figure 1E, parahippocampal cortex, PHC, medial orbitofrontal cortex, mOFC). The straight dynamics also showed a geometric change back and forth along the straight trajectory (Figure 1D, i.e., overlapping continuous line). We note that one neural trajectory was obtained from a population of neurons in this analysis method.
In short, these qualitative observations based on visual inspection suggest that neural population structures may change dynamically through the visual recognition process, and the shift of neural population geometries might occur throughout the cortical and sub-cortical structures across the temporal/frontal/limbic network.
To quantify these geometric patterns occurring at approximately half a second of visual item recognition, we estimated indices for the characteristics of dynamic neural changes in the low-dimensional neural state (Figure 2A). They were as accumulated angle difference weighted by deviance, Σdθ, which reflects a measure similar to the centrifugal force (see STAR Methods for details, Figure 2A, top); mean distance of vectors (d, Figure 2A, top); rotational speed (θ/0.1 s, Figure 2A, bottom); and distance between start and end of trajectory (dS-E, Figure 2A, bottom). We performed bootstrap resampling in each of the 10 neural populations to reconstruct neural geometry as a parental population 1000 times (see STAR Methods).Figure 2Quantitative evaluation of geometric structures according to the rotational features(A) Schematic depictions of the estimation of accumulated angle difference weighted by the deviance, Σdθ. The accumulated angle difference indicates the degree of geometric change in terms of the rotational force across time. Vector distance (d), rotational speed (θ/0.1s), and start to endpoint distance (dS-E) were also estimated.(B) Dendrogram estimated from these four parameter values based on bootstrap resampling across 10 neural populations.(C) Percentage of variance explained by PCA of bootstrap resampling data across 10 neural populations.(D) Clusters detected among the four parameters based on the PCA. Dots represent replicates composed of 20,000 (1000 replicates in each 10 brain regions times two task conditions).(E) Percentage of the identified clusters in each of the 10 brain regions. Each neural population contained two components of neural the best (B) and worst (W) conditions in Exps. 1 and 3, magnitude (M) and delay (D) of the rewards in Exp. 2, and magnitude (M) and probability (P) of rewards in Exp. 4. Colors on the atlas indicate geometry types on visual inspection in Figure 1A.
The bootstrap resampling is an analytic technique to reconstruct the parental population from the observed sample data,^36^ which can provide insight into the types of neural population geometry represented in the parental population. Then, we estimated the parameter values for each replicated neural geometry that provides the parental distribution of the geometric features. This bootstrap resampling was aimed to examine how dominantly a particular neural population geometry was observed in the replicates of 10 neural populations, and thus, we made clustering across all the replicated data to identify whether the single or combination of particular geometric features was observed in some brain regions (otherwise, single cluster occupies all replicates in a single population in a brain region).
We found that across the 10 neural populations, these indices captured geometric features to some extent, corresponding to the rotational geometries observed in the TE and STRt (Figures 2B–2E). For instance, identified clusters based on the dendrogram and principal component analysis (PCA) (Figures 2B and 2C) showed that they possess a rotational characteristic with high rotational speed (Figure 2D, right, red), large Σdθ (Figure 2D, middle, red), large d (Figure 2D, middle, red), and small dS-E (Figure 2D, left, red), which occupy more than 90% of the STRt replicates in the best item condition (Figure 2E, see also Figure 1B, right green trajectory). A smaller rotational structure characterized by smaller values of Σdθ was also captured by another cluster (Figures 2D and 2E, shallow red), which occupied approximately 90% of the STRt replicates in the worst item condition (Figure 2E, see also Figure 1B right, orange trajectory). These rotational features were observed in other temporal brain regions (Figure 2E, see reddish, more than 50% in TE, 20–40% in PRC and PHC), but rarely observed at the frontal/limbic brain regions, such as HPC (less than 10% in all remaining brain regions). In contrast, curvy/straight dynamics were observed in other clusters in the downstream brain regions (Figure 2E, green and blue).
Collectively, across the replicates obtained from the 10 neural populations, a cluster with high rotational speed occupied the STRt and half of the TE populations (Figures 2D and 2E, reddish), while the curvy/straight dynamics occupied most of the replicates in the remaining cortical and subcortical regions (Figures 2D and 2E, blue and green). We note that the meaning of these geometric differences among brain regions will be discussed at the discussion section.
To parameterize these geometric features in more detail, we fitted the Lissajous curve function^37^ to the replicated data, which can mathematically capture all rotational, curvy, and straight geometries by the single function. In the Lissajous function, any two-dimensional geometric features represented by F(x, y) are captured using the following (Equation 1)x=Axcos(ωxt(i)+Φx)+bx(Equation 2)y=Aycos(ωyt(i)+Φy)+bywhere ω and Φ (i.e., ωx, ωy, Φx, and Φy) represent the cycle of the rotation and their deviation as a function of time, t(i). t(i) takes the values from 0 to 0.6 s in all the four experiments, and thus, one cycle of the trajectory is represented as 0 to 3.33 π for ω. For the horizontal and vertical axes, Ax and Ay determine the amplitude of the trajectory, respectively, whereas bx and by determine its location on the PC1-2 plane. In this function, rotational dynamics are represented by the same ω among x and y formulas, and 0.5 π cycle differences in Φ values between x and y formulas (Figure 3A, left). We note that this rotational example represents less than one cycle due to ωx = 3.0 π. In contrast, straight dynamics are represented with the same ω and also same Φ values between the two formulas (Figure 3A, right). Curvy dynamics are represented with some difference of ω and same Φ values (Figure 3A, middle) (see also Figure S4). We fitted this Lissajous curve function to each of the 20,000 bootstrap replicates derived from the 10 neural populations (see STAR Methods, 1000 replicates times 10 populations times two conditions). For instance, three replicated examples obtained from the HPC population were well captured by the Lissajous curve function, as rotational-to-straight trajectories (Figure 3B). We obtained all the estimates for these fitted parameters (Figure 3C), and thereafter, applied clustering to these parameters’ data (Figures 3D–3F) to identify geometry types as a function of the Lissajous curve parameters.Figure 3Quantitative evaluation of geometric structures using the Lissajous curve function(A) Schematic depictions of trajectory geometries using Lissajous function parameters. For all figures, ωx and ωy are 3 π.(B) Three examples of bootstrap replicates for the HPC population fitted by the Lissajous function. These examples were obtained at the best stimulus condition. L indicates the maximum log likelihood. Estimated parameters were as left, ωx, 2.78 π, Φx, −0.11, L, 34.6, ωy, 2.96 π, Φy, 0.10, L, 31.6; middle, ωx, 2.51 π, Φx, −0.03, L, 24.2, ωy, 3.82 π, Φy, −0.08, L, 32.7; right, ωx, 2.21 π, Φx, −0.16, L, 26.9, ωy, 2.78 π, Φy, 0.58, L, 29.1.(C) Probability density estimated for Lissajous parameters obtained from bootstrap replicates across 10 neural populations times two conditions.(D) Dendrogram estimated from Lissajous parameter values based on bootstrap resampling across 10 neural populations times two conditions.(E) Percentage of variance explained by PCA of bootstrap resampling data across 10 neural populations times two conditions.(F) Clusters determined using PCA. Data are shown for PC1 to 3.(G) Reconstructed trajectory in each cluster based on bootstrap resampling. The trajectories in clusters 1–10 were drawn using the median values of the Lissajous parameters in each cluster.(H) Percentage of clusters in each of the 10 brain regions times two conditions. BW: best and worst conditions. MD: magnitude and delay conditions. MP: magnitude and probability conditions. See also Figure S4.
We found that the rotational dynamics (Figures 3G and 3H, clusters C1-C5, reddish) appeared at the TE and STRt, which occupied high percentages of these neural populations (Figure 3H, approximately 70%), and they were also observed in more than 50% of temporal brain regions. Cluster 5 seemed to have the intermediate characteristics between rotational and curvy structures; if we define this cluster 5 as the curvy one, the rotational percentage becomes low in the PRC, PHC, and HPC (30–40%), but not in the STRt (50–70%). Curvy structures were predominantly observed in the CDb population, with more than 40% occurrence (Figures 3G and 3H, Clusters 6 and 7, greenish). These clusters were not observed deterministically but rather stochastically, as also seen in the predominant percentages of intermediate features between curvy and straight dynamics (Figures 3G and 3H, Cluster 8). Straight dynamics were predominant at the frontal brain regions, while they were also observed at the temporal cortices (Figures 3G and 3H, Clusters 8–10, bluish). Even with the STRt replicates, the neural population contained curvy or straight dynamics in more than 20%. These heterogeneous mixtures of replicated signals suggested that neural dynamics emerged in a stochastic manner at population level, with a functional gradient in the temporal/frontal/limbic networks of cortical and subcortical structures. The brain-wide neural population may propagate item information as a heterogeneous mixture at approximately half a second. We note that our results do not reflect trial-by-trial variability of neural population activity. Additionally, we acknowledge that our results contain some noise, which contributes to stochastic differences among replicates to some extent.
Collectively, our meta-analysis results indicated that the rotational, curvy, and straight neural geometries were found in neural populations across temporal-frontal-limbic cortical sub-cortical networks (Figure 1). Orbitofrontal cortex (cOFC and mOFC) and their target subcortical brain region VS predominantly showed curvy, straight, and their intermediate dynamics (Figures 3G and 3H). These observed dynamics exhibited maximum modulation at approximately 0.3 s after visual onset, except in the slowest VS dynamics (Figure 1D). Under our task conditions, the rotational dynamics, in contrast, were observed at the temporal cortices and their connected striatal regions, at the relatively shorter latency around 0.2 s when rotation started (Figure 1B). Whereas different monkeys performed active (Exps. 1 and 3) and passive (Exps. 2 and 4) behavioral tasks, the rotational dynamics were observed predominantly in the STRt and TE (Figure 3H). In contrast, straight dynamics started its geometric change approximately after 0.2 s of the visual onset (Figure 1D, see the distance between initial point S and 0.2 s location), indicating that they follow the rotational dynamics in processing visual information.
In our dataset, monkeys were engaged in the tasks involving stimulus-location association (Exp. 1), stimulus-reward association (Exp. 2), stimulus-reward association for delay and magnitude (Exp. 3), and stimulus-reward association for probability and magnitude (Exp. 4). All the four tasks required monkeys to process visual memory associations to perform the task. Taken together, while these neural geometries are observed during visual recognition of such associations, which reflect visual memory processing, these three geometries were distinctive in terms of geometric patterns and their dynamic changes over time (Figure 4A for summary), in which a rotational/curvy change was followed by a change in straight dynamics.Figure 4Summary of the observed dynamics and anatomical connections in the visual memory pathway(A) Geometries depicted in the same arbitrary scales on the PC1-2 plane for the eight neural populations shown in Figures 1B–1D. The start of the trajectory (S) is aligned to describe each trajectory. e indicates the end of the trajectory at 0.6 s.(B) Proportion of the clusters defined in each of the 10 brain regions are described with the anatomical connection. Reddish: rotational, curvy, straight dynamics. Data from CDh&b and CDb are merged (CD).
Previous studies showed that rotational dynamics have been uncovered broadly in the primary sensory^24^^,^^38^ and motor^39^ cortices, which are closer to the inputs and outputs of the brain, such as motor unit activity.^40^ Other studies have shown that the prefrontal cortex^35^^,^^41^ and parietal cortex^22^^,^^38^ exhibit curvy dynamics. Does this evidence mean that each brain region represents a particular geometry type on the activity of neural population?
It is likely that the brain regions close to the input/output of the central nervous system may represent particular types of geometries, while the higher-order brain regions may represent various types of geometries. Because geometries reflect neural activity change during behavioral tasks, flexibility of the cognitive function must be related to the geometric type-change. In our case, we found curvy dynamics in the CDb, where action information was transferred from the cortices^42^ in the Exp. 3 (Figure 3H), but more heterogeneous geometries were observed from CDh&b in the Exp. 4 (Figure 3H). Although the recorded areas of neurons were not perfectly matched between these two experiments, geometric patterns must be changed related to the task demand to the subjects.
In the present study, we specifically focused on neural dynamics in two core (1) low-dimensional geometries and (2) neural modulation dynamics. First, although these data were obtained from four different laboratories using distinctive behavioral tasks during the passive and active responses of monkeys, the low-dimensional features of neural dynamics are thought to be preserved across mammals in the brain-wide network.^43^^,^^44^ A recent study provides clear evidence that different animal species share and preserve their neural geometries during behavioral tasks.^43^ Indeed, it is suggested that a low-dimensional manifold in the neural state-space might be one of the representational states of biologically relevant information, similar to many combinations of physical properties in the world.^45^ Second, the dynamics in neural modulations examined here are comparable with those using the standard analytic frameworks in the rate coding model, which has provided a huge amount of knowledge corresponding to low-dimensional neural activity modulation in the literature, such as the Gabor function in the visual cortices,^46^^,^^47^ movement direction and muscle force in the motor cortex,^48^^,^^49^ reward value in the parietal cortex^50^ and frontal cortical and subcortical regions,^51^^,^^52^ action values in the striatum,^53^^,^^54^ and comprehension of the location of animals during navigation in the HPC.^55^ Thus, in our analysis, the dynamics of these well-known brain features were compared as the geometric patterns across neural populations during visual recognition.^29^
One concern with our meta-analysis approach is that there may be limitations in data interpretation in terms of data sharing and comparisons across different behavioral tasks and different individual animals. Is it possible to compare the neural population trajectory using accumulated data across animals and tasks with a certain analytical tool? In previous literature, when analyzing neural modulation using a linear regression model, most neurophysiologists compare the extent of neural modulations between different animals and between behavioral tasks. In one study, a challenge was made to compare the neural modulation dynamics as the trajectory geometry between different laboratories’ data.^56^ Thus, the neural trajectory would be comparable among our shared data with greater deliberation.
Our result showed that variability of neural population geometry is distinct among the 10 neural populations examined (Figure 4B, summary for Figure 3H). This raised the question of why such variabilities is observed in the reconstructed neural geometries. Using bootstrap resampling, we reconstructed the parental population of geometry from observed sample data. In this approach, (1) we generated replicates of neural populations through random resampling of neuronal activity, allowing duplicates when resampling from a single neural population dataset (i.e., each population data shown in Table S1). The number of resamples was matched to the number of recorded neurons in each neural population dataset (e.g., 295 and 407 times in TE and PRC, respectively). (2) This process was repeated 1,000 times for each neural population to yield 1,000 replicates. (3) Each replicate was analyzed using dimensionality reduction to produce 1,000 replicated geometries per neural population. This approach provided parental distributions of neural-geometry parameters for each of the 10 neural populations (Figures 2 and 3). Then, clustering analysis revealed that some populations such as, STRt, showed nearly uniform trajectories, while others, such as PRC, displayed multiple trajectory types (Figures 2E and 3H). These trajectory patterns (Figures 3F–3H) indicated that parental populations have considerable variabilities.
The observed variability of these replicated parameters in multidimensional space (Figures 2D and 3F) may reflect functional aspects of neural population dynamics. One possibility is that the variabilities arise from dynamic engagement within the neural population—some neurons contribute to information processing at one moment but may not participate at another moment at the population level (a phenomenon we term “population stochasticity”). Another possibility is that heterogeneity within the neural population contributes to this variability due to the random resampling process. However, the latter explanation is less likely because, with a larger neural population, replicate variability should decrease. Interestingly, such a decrease did not occur even in the largest sampled population, HPC, with 590 neurons (Figure 3, HPC).
These findings suggest that neural population stochasticity is key to understanding how large-scale networks process information. Importantly, we note that this variability does not reflect trial-by-trial variability of neural population activity, as no simultaneous neural recordings were made in this study.
Our findings would add to the emerging literature describing how visual inputs alter the brain-wide neural dynamics associated with visual memories by connecting neural geometry types and their alignments across many brain regions. Previous studies have shown the existence of different types of neural population dynamics in each individual study.^22^^,^^24^^,^^35^^,^^38^^,^^39^^,^^41^ Although some of these dynamics may reflect task demand, as observed in the dorsolateral prefrontal cortex,^35^^,^^57^ they are difficult to disentangle from changes in behavioral and neural activity levels, and may involve some transformation of information for behavioral responses.^58^ It is possible that the dorsal and motor-related brain regions have this type of flexibility in their dynamics, as partly observed in this study in the CDb, where curvy dynamics were predominantly observed (Figure 4B).
Our results raise the possibility that geometric features determine the important neural mechanisms widely observed in the brain. For instance, the stochastic gradient relative distance to the visual input may reflect the dynamics of the neural circuitry in half a second (Figure 4B). In contrast, STRt deterministically showed rotational geometries in the bootstrap replicates (Figure 4B, reddish clusters occupy almost 70% of replicates). The difference between stochastic/deterministic observations of geometry indicates two possible functional meanings. One possibility is that neural population may process information with curvy dynamics in some trials, but in other trials they process information with rotational dynamics. The other possibility is that homogeneity/heterogeneity of neuronal activity in a population may contribute to the variability in the replicates, though this might not be likely because bootstrap resampling must provide more similar replicates if the number of observations is large (see, Hippocampal results, Table S1, n = 590 neurons which is the largest population).
For summary of the present study (Figure 4B), the unidimensional straight dynamics in the hippocampal-frontal circuitry may reflect memory access during visual recognition, such as location and reward. The rotational dynamics in the TE/STRt might reflect the visual recognition process, during which recurrent feedback signals change the circuit dynamics. Future studies should test the underlying functional mechanisms and define whether the engagement of the change in behavior and/or task context is best considered for whole-brain neural population activity. Regardless of the mechanism, the shift of modulation structures in the lower-dimensional neural space could play a fundamental role in brain-wide information processing, such as transforming visual feature recognition to memory access.
Our result was meta-analysis, and thus, cannot be applied to understand individual cases.
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Hiroshi Yamada (h-yamada@md.tsukuba.ac.jp).
No materials available in this study.
All data and analysis codes used in this study are available in the supporting files.
The authors thank Takashi Kawai, Yoshiko Yabana, Yuki Suwa, and Shiho Nishino for their technical assistance. We appreciate Shigeru Shinomoto for his comments. Monkey FU was provided by NBRP “Japanese Monkeys” through the National Bio Resource Project of the 10.13039/501100001700MEXT, Japan. Funding: This research was supported by 10.13039/501100001691JSPS KAKENHI (Grant Numbers JP:15H05374, 22H04832), JST Moonshot R&D JPMJMS2294 (H.Y.), and the 10.13039/501100001809National Natural Science Foundation of China (Grant 32271088) (Y.N.).
H.Y. conceptualized the study. H.Y., Y.N., T.M., O.H., and J.K. designed the experiments. H.Y., H.C., J.K., Y.H., Y.I., and T.M. performed the experiments. H.Y. and Y.T. developed the analytical tools. H.Y., H.C., J.K., and Y.H. analyzed the data. H.Y., H.C., J.K., Y.H., and Y.T. wrote the manuscript. All authors edited and approved the final version of the manuscript.
The authors declare no competing of interests.
REAGENT or RESOURCESOURCEIDENTIFIERSoftware and algorithmsR software 4.4R projecthttps://www.r-project.org/MATLAB 2020bMathWorks Inc.https://mathworks.com/products/matlab.htmlAdobe Illustrator CS6Adobehttps://www.adobe.com/products/illustrator.htm
Subjects and experimental procedures are as follows. Nine rhesus monkeys were used in the present study (Exp. Macaca mulatta, A, 9.3 kg, male; Macaca mulatta, D, 9.5 kg, male; Exp Macaca mulatta, WK, 12.0 kg, male; Macaca mulatta, SP, 7.0 kg, male; Exp Macaca mulatta, BI, 8.2 kg, male; Macaca mulatta, FG, 11.0 kg, male; Macaca mulatta, ST, 5.2 kg, male; Exp Macaca mulatta, SUN, 7.1 kg, male; Macaca fuscata, FU, 6.7 kg, female). All experimental procedures were approved by the Institutional Animal Care and Use Committee of Laboratory Animals approved by Peking University (Exp. 1, project number Psych-YujiNaya-1), and Animal Care and Use Committee of the National Eye Institute, and complied with the Public Health Service Policy on the Humane Care and Use of Laboratory Animals (Exp. 2, protocol number NEI-622), the Animal Ethics Committee of the National Institutes for Quantum Science and Technology (Exp. 3, protocol no. 11-1038-11), and the Animal Care and Use Committee of the University of Tsukuba (Exp. 4, protocol no 23-057). All procedures were performed in reference to the US Public Health Service Guide for the Care and Use of Laboratory Animals.
The animals performed the task under dim light conditions in an electromagnetically shielded room. The task started with an encoding phase, which was initiated by the animal pulling a lever and fixating on a white square (0.6°) presented within one of four quadrants at 12.5° (monkey A) or 10° (monkey D) from the center of the touchscreen (3M™ MicroTouch™ Display M1700SS, 17 in), situated approximately 28 cm from the subjects. The eye position was monitored using an infrared digital camera with a sampling frequency of 120 Hz (ETL-200, ISCAN). After fixation for 0.6 s, one of six items (3.0° for monkey A and 2.5° for monkey D, radius) was presented in the same quadrant as a sample stimulus for 0.3 s, followed by another 0.7 s fixation on the white square. If the fixation was successfully maintained (typically < 2.5°), the encoding phase ended with the presentation of a single drop of water.
The encoding phase was followed by a blank interphase delay interval of 0.7–1.4 s during which no fixation was required. The response phase was initiated using a fixation dot presented at the center of the screen. One of the six items was then presented at the center for 0.3 s, as a cue stimulus. After another 0.5 s delay period, five disks were presented as choices, including a blue disk in each quadrant and a green disk in the center. When the cue stimulus was the same as the sample stimulus, the animal was required to make a choice by touching the blue disk in the same quadrant as the sample (i.e., the match condition). Otherwise, the subject was required to choose the green disk (i.e., non-match condition). If the animal made the correct choice, four–eight drops of water were provided as a reward; otherwise, an additional 4 s was added to the standard inter-trial interval (1.5–3 s). The number of reward drops was increased to encourage the animal to maintain good performance in the latter phase of a daily recording session, which was typically conducted in blocks (e.g., a minimal set of 60 trials with equal numbers of visual items presented in a match/nonmatch condition). During the trial, a large gray square (48° on each side) was presented at the center of the display as a background. At the end of the trial, all stimuli disappeared, and the entire screen displayed a light red color during the inter-trial interval. The start of a new trial was indicated by the reappearance of a large gray square on the display, at which point the monkey could pull the lever, triggering the appearance of a white fixation dot.
In the match condition, sample stimuli were chosen pseudo-randomly from six well-learned visual items, and each item was presented pseudo-randomly within four quadrants, resulting in 24 (6 × 4) configuration patterns. In the non-match condition, the location of the sample stimulus was randomly chosen from the four quadrants, and the cue stimulus was randomly chosen from the remaining five items that differed from the sample. The match and non-match conditions were randomly presented in a ratio of 1, resulting in 30 (24 + 6) configuration patterns. The same six stimuli were used during all the recording sessions.
Animals learned the scene-object associations. After the monkeys fixated on the red-square fixation point on the scene image for 0.6–1 s, the fixation cue disappeared, and two visual items (objects of different values) appeared simultaneously in a different hemifield (for training and neuronal testing) or the same hemifield (for pharmacological experiments). A reward was given after the monkeys made a saccade to the stimulus and maintained fixation for 0.2 s. Half of the fractal visual items were associated with a large reward (0.3 mL), and the other half were associated with a small reward (0.1 mL). This reward association changed depending on the scene (Figure S1D). Passive Viewing Task. One of the two scene images was presented for 0.8 s randomly. If the monkey fixated on a central red square, two to four fractals were presented sequentially on the scene image within the neuron’s receptive field (presentation time, 0.4 s; interstimulus interval, 0.4 s; Figure S1C). A liquid reward (0.2 mL) was delivered 0.3 s after the last object was presented. Thus, reward occurrence was not associated with any of the visual items. Each item was presented at least seven times per session.
The monkeys were seated on a primate chair inside a dark, sound-attenuated, electrically shielded room. A touch-sensitive bar was mounted on the chair. The visual stimuli were displayed on a computer video monitor placed in front of the animals. Each of the six cues was associated with a combination of reward size (1 drop; 3 or 4 drops) and reward delay (0, 3.3, and 6.9 s). The trials began when the monkey touched the bar. A visual cue appeared, and the monkey released a bar when a red spot (waiting signal) turned green (go signal) after a variable interval. If the monkey released the bar 0.2–1 s after this go signal, the trial was considered correct and the spot turned blue (correct signal). A liquid reward of a small (1 drop, approximately 0.1 mL) or large amount (3 drops, except for monkey BI, 4 drops) was delivered immediately (0.3 ± 0.1 s) or with an additional delay of either 3.3 ± 0.6 s or 6.9 ± 1.2 s after correct release of the bar. The cues were chosen with equal probability and were independent of the preceding reward condition. Anticipatory bar releases (before or no later than 0.2 s after the appearance of the go signal) and failure to release the bar within 1 s of the appearance of the go signal were counted as errors. In the error trials, the trial was terminated immediately, all visual stimuli disappeared, and following inter trial interval (1 s), the trial was repeated; that is, the reward size/delay combination remained the same as that in the error trial. Behavioral control and data acquisition were performed using a real-time experimentation system (REX).^59^ The Neurobehavioral Systems Presentation software was used to display the visual stimuli (Neurobehavioral Systems).
The animals performed one of two visually cued lottery a single-cue or a choice task. Neuronal activity was only recorded during the single-cue task.
Animals performed the task under dim lighting conditions in an electromagnetically shielded room. Eye movements were measured using a video camera system at 120 Hz (EyeLink, SR Research). Visual stimuli were generated using a liquid-crystal display at 60 Hz, placed 38 cm from the monkey’s face when seated. At the beginning of the single-cue task trials, the monkeys had 2 s to align their gaze within 3° of a 1°-diameter gray central fixation target. After a fixation for 1 s, a pie chart was presented for 2.5 s, to provide information regarding the probability and magnitude of rewards in the same location as the central fixation target. The probability and magnitude of the rewards were associated with the number of blue and green 8° pie chart segments, ranging from 0.1 to 1.0 mL in 0.1 mL increments for magnitude, and 0.1 to 1.0 in 0.1 increments for probability. Following a 0.2 s interval from the removal of the pie chart, a 1 kHz or 0.1 kHz tone of 0.15 s duration was provided to indicate reward or no-reward outcomes, respectively. After a 0.2 s interval following the high tone, a fluid reward was delivered, whereas no rewards were delivered following the low tone. An inter-trial interval of 4–6 s was used. During the choice task, animals were instructed to choose one of two peripheral pie charts, each of which indicated either the probability or magnitude of an upcoming reward. The two target options were presented for 2.5 s at 8° to the left or right of the central fixation location. The animals received a fluid reward as indicated by the green pie chart of the chosen target, with the probability indicated by the blue pie chart. Otherwise, no reward was delivered.
A total of 100 pie charts composed of 10 levels of probability and magnitude of rewards were used in the experiments. In the single-cue task, 100 pie charts were presented once in random order. In the choice task, two pie charts were randomly assigned to the two options. During one electrophysiological recording session, approximately 30–60 trial blocks of the choice task were interleaved with 100–120 trial blocks of the single-cue task.
To record the single-unit activity, we used a 16-channel vector array microprobe (V1 X 16-Edge, NeuroNexus), 16-channel U-Probe (Plexon), tungsten tetrode probe (Thomas RECORDING), or single-wire tungsten microelectrode (Alpha Omega). Electrophysiological signals were amplified, bandpass-filtered (200–6000 Hz), and monitored. Single-neuron activity was isolated based on spike waveforms, either online or offline. For both clustering and offline sorting, the activities of all single neurons were sampled when the activity of an isolated neuron demonstrated a good signal-to-noise ratio (>2.5). The signal-to-noise ratio was visually checked by calculating the range of background noise against the spike amplitude, which was monitored online using the OmniPlex Neural Data Acquisition System, or offline using the sorter software Plexon. The recorded neurons were not blinded. The sample sizes required to detect the effect sizes (numbers of recorded neurons, recorded trials in a single neuron, and monkeys) were estimated based on previous studies.^32^^,^^60^ Neural activity was recorded during 60–240 trials of the ILR task. We recorded 590 hippocampal neurons, among which the recording sites appeared to cover all subdivisions (i.e., the dentate gyrus, CA3, CA1, and subicular complex).
We used conventional techniques to record the single-neuron activity in the STRt, including the caudate and putamen tails. A tungsten microelectrode (1–3 MΩ FrHC; 0.5-1.5 MΩ Alpha Omega Engineering) was used to record single-neuron activity. The recording site was determined using a grid system that allowed electrode penetration at 1 mm intervals. We amplified and filtered (0.3 to 10 kHz; Model 1800, A-M Systems; Model MDA-4I, BAK) signals obtained from the electrodes and collected at 1 kHz. Single neurons were isolated online using custom voltage–time window discriminator software (Blip; available at http://www.robilis.com/blip/). The presumed medium spiny neurons were identified based on their low baseline activity (<3 spikes/s) and broad action potentials.^61^ The recorded neurons were not blinded. The sample sizes required to detect the effect sizes (numbers of recorded neurons, recorded trials in a single neuron, and monkeys) were estimated based on previous studies.^62^^,^^63^ Neural activity was recorded during 10–30 trials of the passive viewing task. We recorded 115 medium spiny neurons in the STRt. In Exp. 2, only a single-neuron recording was performed online. We note that we termed the scene and object for two visual stimuli in our previous study,^64^ but here we termed them scene and item.
Conventional techniques were used to record single-neuron activity in the dorsal part of the head of the caudate nucleus (CD). A tungsten microelectrode (1.1–1.5 MΩ, Microprobes for Life Science; 1.0 MΩ, Alpha Omega Engineering Ltd.) was used to record single-neuron activity. The electrophysiological signals were amplified and monitored using a TDT recording system (RZ2, Tucker-Davis Technologies, TDT). Single-neuron activity was manually isolated based on the online spike waveforms. The activity of all single neurons was sampled from the activity of presumed projection neurons, which are characterized as having a low spontaneous discharge rate (<2 spikes/s) outside the task context and exhibiting phasic discharges in relation to one or more behavioral task events.^61^ Neural activity was recorded during 100–120 trials per block in the delayed-reward task. We recorded the CD of the left or right hemisphere in each of the three monkeys in the experiment, with 150 CD neurons (51, 31, and 68 from the BI, FG, and ST, respectively).
Conventional techniques were used to record single-neuron activity in the DS, VS, cOFC (area 13M), and mOFC (area 14o). A tungsten microelectrode (1–3 MΩ, FHC) was used to record single-neuron activity. Electrophysiological signals were amplified, band-pass filtered (50–3,000 Hz), and monitored using a TDT recording system (RZ5D, Tucker-Davis Technologies, TDT). Single-neuron activity was manually isolated based on the online spike waveforms. The activity of all single neurons was sampled when the activity of an isolated neuron demonstrated a good signal-to-noise ratio (>2.5). The signal-to-noise ratio was calculated online as the ratio of the spike amplitude to the baseline voltage range on the oscilloscope. The recorded neurons were not blinded. The sample sizes required to detect the effect sizes (numbers of recorded neurons, recorded trials in a single neuron, and monkeys) were estimated based on previous studies.^51^^,^^53^^,^^65^ Neural activity was recorded during 100–120 trials of the single-cue task. Neural activity was not recorded during selection trials. We recorded the neurons of a single right hemisphere in each of the two 194 DS neurons (98 and 96 from monkeys SUN and FU, respectively), 144 VS neurons (89 SUN and 55 FU), 190 cOFC neurons (98 SUN and 92 FU), and 158 mOFC neurons (64 SUN and 94 FU). In Exp. 1, only a single-neuron recording was performed online. We recorded presumed medium spiny projection neurons from the DS and VS^23^^,^^54^^,^^66^^,^^67^^,^^68^ and presumed pyramidal neurons in the cOFC and mOFC.^23^^,^^51^^,^^52^
For statistical analysis, we used the statistical software package MATLAB (MathWorks, Exps. 1 and 2), and R (Exps. 3 and 4) for conventional analyses such as linear regression and ANOVA. To analyze the regression matrix using PCA, we used R software. All statistical tests for the neural analyses were two-tailed.
No new behavioral results were included; however, the procedure for the behavioral analysis was as follows.
We previously reported that two monkeys learned to retain the item and location information of a sample stimulus.^69^ Here, we describe the analysis steps used to check whether the monkey used both item and location information to perform the task.
To examine this, we compared the animals' actual correct rates during the recording to random correct rates (chi-square test). The ILR response phase had five options, resulting in a 20% random correct rate. If the animal used an incorrect strategy, such as only retaining the location information of the sample stimulus and ignoring the item information, the correct rate for the match condition would be 100% and that for the nonmatch condition would be 0. Based on the above considerations, we examined the correct rates of the two animals in the match and nonmatch conditions, respectively. In general, the average correct rates for both animals in the match and nonmatch conditions were well above chance levels after training.
We previously reported that two monkeys switched their behavior depending on the value of the item based on the scene.^64^ Here, we describe how to check whether the monkey learned both the scene and item information. We calculated the correct rate for the scene-based object-value task. Because the two scenes appear in random sequences, the monkey must switch object choice if the scene has changed. After performing more than 160 trials, the correct rate reached a plateau above chance. The monkey was able to switch object choices immediately after the scene changed. Once the monkeys learned this extensively, their choice behavior became automatic, as the choice tended to occur even when the reward was not delivered after saccades to high-valued items according to the scene.
We previously reported that the three monkeys behaved based on temporally discounted values that integrated both delay and reward size information provided by visual stimuli.^70^ Here, we describe an analysis to check how monkeys discount reward values for delay and reward information. Error rates in task performance were calculated by dividing the total number of errors by the total number of trials for each reward condition and then averaged across all sessions. The average error rates were fitted to the inverse function of reward size with hyperbolic temporal E = 1+kD/ aR (E: average error rates, D: delay, R: reward size, k: discounting factor, a: incentive impact), and exponential temporal E = e^-kd^/aR. We used the ‘optim’ function in R, evaluated the goodness of fit of the two models by least-squares minimization, and compared the models by leave-one-out cross-validation as described previously (Minamimoto et al.,^71^ 2009).
We previously reported that monkey behavior depends on expected values, defined as the probability time magnitude.^23^ We described the analysis steps to check whether the monkey’s behavior reflected task parameters, such as reward probability and magnitude. Importantly, we showed that the monkey’s choice behavior reflected the expected values of the rewards, that is, the probability multiplied by the magnitude. For this purpose, the percentage choosing the right-side option was analyzed in the pooled data using a general linear model with a binomial (Equation 3)PchoosesR=1/(1+e−z)where the relationship between PchoosesR and Z is given by the logistic function in each of the following three number of pie segments (M1), probability and magnitude (M2), and expected values (M3).(Equation 4)M1:Z=b0+b1NpieL+b2NpieRwhere b0 is the intercept, and NpieL and NpieR are the number of pie segments contained in the left and right pie chart stimuli, respectively. The values of b0 to b2 are free parameters and estimated by maximizing the log likelihood.(Equation 5)M2:Z=b0+b1PL+b2PR+b3ML+b4MRwhere b0 is the intercept; PL and PR are the probabilities of rewards for the left and right pie chart stimuli, respectively; and ML and MR are the magnitudes of rewards for the left and right pie chart stimuli, respectively. The values of b0 to b4 are free parameters and estimated by maximizing the log likelihood.(Equation 6)M3:Z=b0+b1EVL+b2EVRwhere b0 is the intercept and EVL and EVR are the expected values of rewards as probability multiplied by magnitude for the left and right pie chart stimuli, respectively. The values of b0 to b2 are free parameters and estimated by maximizing the log likelihood. We identified the best model to describe monkey behavior by comparing goodness-of-fit based on Akaike’s information criterion and Bayesian information criterion.^72^
Peri-stimulus time histograms were constructed for each single-neuron activity aligned at the onset of the visual stimulus. Average activity curves were smoothed for visual inspection using a Gaussian kernel (σ = 20, 15, 10, and 50 ms in Exps. 1–4, respectively), whereas the Gaussian kernel was not used for statistical tests.
To ensure that the four different datasets were as fair as possible, we used the same criteria to analyze the neural activity. For the neural analyses, we used the following four 1) the same analysis window size, 2) visual response within a short time (0.6 s), 3) neural modulations detected at the same significance level (P < 0.05), and 4) a general linear model (ANOVA in Exps. 1 and 2 and the linear regression in Exp. 3 and 4). The details of these analytical procedures for the rate coding and dynamic models are shown below.
For neural responses during the encoding phase after the sample presentation, we evaluated the effects of “item” and “location” for each neuron using two-way ANOVA (P < 0.05 for each). We analyzed neurons that were tested in at least 60 trials (10 trials for each stimulus and 15 trials for each location). On average, we tested 100 trials for each neuron. These results have been previously reported.^69^
For neural responses during the appearance of the visual item, we evaluated the effects of “item” and “scene” for each neuron using paired t-test (P < 0.05 with Bonferroni correction). These results have been previously reported.^64^
The neural discharge rates (F) were fitted using a linear combination of the following (Equation 7)F=b0+bdDelay+bmMagnitudewhere Delay and Magnitude are the delay and magnitude of the reward, respectively, as indicated by the visual stimulus. b0 is the intercept. If bd and bm were not zero at P < 0.05, the discharge rates were regarded as being significantly modulated by that variable. These results have been previously reported.^70^
The neural discharge rates (F) were fitted using a linear combination of the following (Equation 8)F=b0+bpProbability+bmMagnitudewhere Probability and Magnitude are the probability and magnitude of the rewards, respectively, as indicated by the pie chart. b0 is the intercept. If bp and bm were not zero at P < 0.05, the discharge rates were regarded as being significantly modulated by that variable. These results have been previously reported.^23^
We analyzed neural activity during an identical 0.6 s duration from the sample onset (Exp. 1), item onset (Exp. 2), CUE onset (Exp. 3), and cue onset (Exp. 4). To obtain a time series of neural firing rates within this time period, we estimated the firing rates of each neuron for every 0.05 s time bin (without overlap) during the analysis periods. A Gaussian kernel was not used.
We used a general linear model to determine how items and locations (Exp. 1), items and scenes (Exp. 2), delay and magnitude of rewards (Exp. 3), and the probability and magnitude of the rewards (Exp. 4) affect the activity of each neuron in the neural populations. Each neural population was composed of all the recorded neurons in each brain region.
First, we set six visual items and four locations as categorical variables. We then described the average firing rate of neuron i at time t as a linear combination of the item and the location in each neural (Equation 9)F(i,t,k)=b0(i,t)+b1(i,t)Item(k)+b2(i,t)Location(k)where F(i,t,k) is the average firing rate of neuron i at time t in trial k, Item~(k)~ is the types of items cued to the monkey in trial k, and Location~(k)~
is the types of locations cued to the monkey in trial k. The regression coefficients b0(i,t), b1(i,t), and b2(i,t) describe the degree to which the firing rates of neuron i depend on the mean firing rates (hence, firing rates independent of task variables, item, and location), the degree of firing rate in each item relative to the mean firing rates, and the degree of firing in each location relative to the mean firing rates, respectively, at a given time t during the trials. The interaction term is not included in the model.
In the analysis, we performed preference ordering for item and location in each neuron. Item~(k)~ and Location~(k)~ were rank-ordered items and locations, respectively, cued to the monkey in trial k. Items 1–6 and locations 1–4 were rank-ordered from the most preferred to least preferred, respectively, defined as the mean firing rate during the entire analysis time window from 0.08 to 0.6 s. This preference ordering did not change over time t for each neuron n.
We first set eight items and two scenes as the categorical variables. We then described the average firing rate of neuron i at time t as a linear combination of the item and scene in each neural (Equation 10)F(i,t,k)=b0(i,t)+b1(i,t)Item(k)+b2(i,t)Scene(k)where F(i,t,k) is the average firing rate of neuron i at time t in trial k, Item~(k)~ is the types of items cued to the monkey in trial k, and Scene~(k)~ is the types of scene stimuli cued to the monkey in trial k. The regression coefficients b0(i,t), b1(i,t) and b2(i,t) describe the degree to which the firing rates of neuron i depend on the mean firing rates (hence, firing rates independent of task variables, item and scene), the degree of firing rate in each item relative to the mean firing rates, and the degree of firing in each scene relative to the mean firing rates, respectively, at a given time t during the trials. The interaction term was not included in the model.
In the analysis, Item~(k)~ and Scene~(k)~ were the rank-ordered item and scene, respectively, cued to the monkey in trial k. Items 1 to 8 and Scenes 1 and 2 were rank-ordered from the most preferred to least preferred, respectively, defined as the mean firing rate during the whole analysis 0.6 s window after the item onset. This preference ordering did not change over time t for each neuron n.
We first set the delay and magnitude as 0, 3.3, and 6.9 s and one and three drops of rewards, respectively, during the behavioral task. In the analysis, we normalized these values from 0 to 1 divided by the maximum values in 0, 0.48, and 1 for delay, and 0.33, 0.66, and 1 for magnitude. This is because these values affect the extent of the regression subspace between two continuous variables. We then described the average firing rate of neuron i at time t as a linear combination of the delay and magnitude in each neural (Equation 11)F(i,t,k)=b0(i,t)+b1(i,t)Delay(k)+b2(i,t)Magnitude(k)where F(i,t,k) is the average firing rate of neuron i at time t in trial k, Delay~(k)~ is the normalized delay to obtain a reward cued to the monkey in trial k, and Magnitude (k)
is the normalized number of reward drops cued to the monkey in trial k. The regression coefficients b0(i,t) to b2(i,t) describe the degree to which the firing rates of neuron i depend on the mean firing rates (hence, firing rates independent of task variables), delay in rewards, and magnitude of rewards, respectively, at a given time t during the trials.
We first set the probability and magnitude as 0.1 to 1.0 and 0.1 to 1.0 mL, respectively. We did not normalize these values because they were originally prepared from 0 to 1. We then describe the average firing rate of neuron i at time t as a linear combination of probability and magnitude in each neural (Equation 12)F(i,t,k)=b0(i,t)+b1(i,t)Probability(k)+b2(i,t)Magnitude(k)where F(i,t,k) is the average firing rate of neuron i at time t in trial k, Probability~(k)~ is the probability of the reward cued to the monkey in trial k, and Magnitude~(k)~
is the magnitude of the reward cued to the monkey in trial k. The regression coefficients b0(i,t) to b2(i,t) describe the degree to which the firing rates of neuron i depend on the mean firing rates (i.e., firing rates independent of task variables), probability of rewards, and magnitude of rewards, respectively, at a given time t during the trials.
We used the regression coefficients (i.e., the regression table in the case of ANOVA) described in Equations 9, 10, 11, and 12 to identify how the dimensions of the neural population signals were composed of information related to the item and location (Exp. 1), item and scene (Exp. 2), delay and magnitude (Exp. 3), and probability and magnitude (Exp. 4) as aggregated properties of individual neural activity. In this step, an encoding model is constructed in which the regression coefficients are explained by a temporal structure in the neural modulation of two categorical variables (Exps. 1 and 2), or two continuous variables (Exps. 3 and 4) at the population level. Our procedures involve targeted dimensionality reduction using the regression subspace^19^ and are aimed at describing neural modulation dynamics.^29^
We used PCA to identify the dimensions of the neural population signal in orthogonal spaces composed of two variables in each neural population of the four experiments. For each neural population, we first prepared a two-dimensional data matrix X of size N(n)×M
(C×T). The regression coefficient vectors b1(i,t) and b2(i,t) in Equations 9, 10, 11, and 12, whose rows correspond to the total number of neurons (n) in each neural population and columns correspond to C, the total number of conditions (that is, six items and four locations in Exp. 1, eight items and two scenes in Exp. 2, delay and magnitude in Exp. 3, and probability and magnitude in Exp. 4), and T is the total number of analysis windows (i.e., 0.6 s divided by the window size bin, 0.05 s, 12 bin). A series of eigenvectors was obtained by applying PCA once to the data matrix X in each neural population. The PCs of this data matrix are vectors v(a) of length N(n) and the total number of recorded neurons if M
(C×T) > N(n); otherwise, the length is M
(C×T). PCs were indexed from the principal components and explained the most to least variance. The eigenvectors were obtained using the prcomp () function in R software. We did not include the intercept term b0(i,t) to focus on the neural modulation by the variables of interest.
When we applied PCA to data matrix X, we decomposed the matrix into eigenvectors and eigenvalues. Each eigenvector had a corresponding eigenvalue. In our analysis, the eigenvectors at time t represented a vector, for example, in the space of delay and magnitude in Exp. 3. The eigenvalues at time t for the delay and magnitude were scalars, indicating the extent of variance in the data in that vector. Thus, the first PC was the eigenvector with the highest eigenvalue. We analyzed the eigenvectors for the top two PCs (PC1 and PC2) in the following analyses to describe the geometry in the most predominant dimension. PCA was applied once to each neural population; thus, the total variance contained in the data differed among the neural populations.
To examine the significance of the population structures described by PCA, we performed three shuffle controls. The two-dimensional data matrix X was randomized by shuffling in three ways. In shuffled control 1, matrix X was shuffled by permutating the allocation of neuron n at time i. This shuffle provided a data matrix X of size N(n)×M
(C×T), eliminating the temporal structure of neural modulation by condition C in each neuron but retaining the neural modulations at time t at the population level. In shuffled control 2, matrix X was shuffled by permutating the allocation of time i in each neuron n. This shuffle provided a data matrix X of size N(neuron)×M
(C×T), eliminating the neural modulation structure under condition C maintained in each neuron but retaining the neural modulation in each neuron at the population level. In shuffled control 3, matrix X was shuffled by permutating the allocation of both time i and neuron n. In these three shuffle controls, matrix X was estimated to be 1,000 times. PCA performance was evaluated by constructing the distributions of the explained variances for PC1 to PC12. The statistical significance of the variances explained by PC1 and PC2 was estimated based on the 95th percentile of the reconstructed distributions of the explained variance or bootstrap standard errors (i.e., standard deviation of the reconstructed distribution). We note that because the significant dimensions of neural populations dynamics differed the 10 neural populations, we analyzed the neural dynamics at the top two dimension, PC1 and 2.
We detected roughly three different types of neural rotational, curvy, and straight dynamics (Figure S4), as interpreted based on the patterns observed in our analysis and inspired by concepts in prior studies of neural dynamics. Rotational dynamics was like circle, and curvy one was like half cycle, and straight dynamics was like line. We did not include stable dynamics, in which trajectory geometry stay a fixed point.
We evaluated the characteristics of the eigenvectors for PC1 and PC2 (i.e., geometry) in each neural population in terms of vector angle, size, and deviation. The eigenvectors were evaluated for each of the task parameters described item and location in Exp. 1, item and scene in Exp. 2, delay and magnitude in Exp. 3, and probability and magnitude in Exp. 4. The angle is the vector angle from the horizontal axis from 0° to 360° against the main PCs. The size is the length of the eigenvector. The deviation is the difference between the vectors (i.e., distance between two vectors). The deviation from the mean vector for each neural population was estimated. These three eigenvector characteristics were compared among the populations at P < 0.05, using the Kruskal–Wallis test and Wilcoxon rank-sum test with Bonferroni correction for multiple comparisons as the basic analyses.^29^ The vector during the first 0.1 s was extracted from these basic analyses.
To evaluate the trajectory geometry using their selected feature, we estimated the accumulated angle difference weighted by the (Equation 13)∑t=St=Edθwhere the d is deviation between the vectors at times t and t+1, θ is the angle difference between vectors at times t and t+1, S is zero, and E is the time period to stop the estimation, i.e., 0.6 s. This index is analogous to the rotational force accumulated over time. If the value of the accumulated angle difference was close to zero, the population geometry was stable, such as a straight or non-dynamic structure, that is, it remained at some point in the PC1-2 plane. In addition to Σdθ, d (distance between vectors), rotational speed Σθ/0.1s, and ds-e, such as start to end distance were estimated.
To evaluate the trajectory geometry without the selected features, we used the Lissajous curve function, which describes any geometric pattern in a plane using F(x,y):(Equation 14)x=Axcos(ωxt(i)+Φx)+bx(Equation 15)y=Aycos(ωyt(i)+Φy)+bywhere ω and Φ represent cycle of the rotation and their deviance as a function of time, t(i). Ax and Ay represent the amplitudes of the trajectory, whereas bx and by represent the location of the trajectory. For ω, 3.33 π indicates that one cycle since the analysis window is 0.6 s. Φ is 0 to 2 π for one cycle. We estimated ωx, Φx, bx, ωy, and Φy, by parameters by estimating maximum loglikelihood of the model. Nonlinear least squares in the nls() function in the R program was used. A time series of eigenvectors for PC1 and PC2 in a 0.05 s analysis windows (12 data points) were used with a sliding average between three time points (hence, 0.15 s time resolution).
We estimated Σdθ, mean d, rotational speed Σθ/0.1s, and ds-e, such as start to end distance using a parametric bootstrap resampling method.^36^ In each neural population, the neurons were randomly resampled in duplicate, and a data matrix X of size N(neuron)×M
(C×T) was obtained. PCA was applied to the data matrix X. Note that the X is regression matrix. The time series of eigenvectors was obtained, and these four features were estimated from the neural trajectory. This resampling was conducted 1,000 times in each neural population, and the distributions of these four parameters were obtained.
Following the bootstrap resampling, we applied clustering of these parameters based on PCA and a dendrogram across the replicates in the 10 brain regions, such as 20,000 replicates (10 brain regions times two conditions times 1,000 replicates). Based on this clustering, proportion of the identified clusters in each brain region was estimated. We note that bootstrap resampling and clustering across all the 10 brain region’s replicates allow us to identify how strongly a particular neural population geometry was observed in each of the 10 neural populations.
The Lissajous curve parameters for the replicated trajectory were estimated using a bootstrap resampling method.^36^ In each neural population, the neurons were randomly resampled in duplicate, and a data matrix X of size N(neuron)×M
(C×T) was obtained. Note that the X is regression matrix. PCA was applied to the data matrix X. The time series of eigenvectors were obtained for PC1 and PC2, which describe the trajectory. The fitted parameters using the Lissajous curve function were estimated using the nls() function in R program. This resampling was conducted 1,000 times in each neural population, and the distributions of the Lissajous parameters were obtained.
Following the bootstrap resampling, we applied clustering of these parameters based on PCA and a dendrogram across the replicates in the 10 brain regions, such as 20,000 replicates (10 brain regions times two conditions times 1,000 replicates). In this process, the omega ratio (ωx/ωy) and phi difference (Φx-Φy) were also used, in addition to the ωx, ωy, Φx, and Φy. Based on this clustering, proportion of the identified clusters in each brain region was estimated. We used the median of the estimated parameters in a cluster to describe the trajectory geometries. We note that bootstrap resampling and clustering across all the 10 brain region’s replicates allow us to identify how strongly a particular neural population geometry was observed in each of the 10 neural populations.
Unless otherwise indicated, all data were presented as means with a distribution (Boxplot or data plot). The statistical analyses performed were indicated in the main text and detailed in STAR Methods. The statistical comparisons were made in R software, including shuffle control and bootstrap resampling.