Behavioral state coding by molecularly defined paraventricular hypothalamic cell type ensembles

1Neuron activity and neuronal gene expression are considered distinct facets of brain function. Large-scale neuron activity recordings reveal discrete groups of neurons that respond similarly, and gene expression brain atlases show differential expression patterns across the brain. Single-cell transcriptomic studies reveal that the brain is composed of hundreds of molecularly defined neuronal subtypes. Molecular markers facilitate repeatedly returning to the same type of neurons and also serve as genetic elements that can be used for targeting functionally similar neurons for causality testing or potentially therapeutics. These applications can only be justified if molecular markers categorize neurons into groups with similar neuronal dynamics that together encode distinct behavioral states. However, there is limited knowledge about the relationship between activity patterns of cells in the neural ensemble and their underlying gene expression profiles. Therefore, the merit of using neuronal gene expression as a proxy for neurons with similar activity patterns is controversial.

2We set out to determine (i) the role of molecularly defined cell types for encoding different behavioral states and (ii) the association of neuronal marker genes with distinct functional response types of the neuronal ensemble from the mouse paraventricular hypothalamus (PVH). The PVH mediates appetite, stress, thirst, and autonomic functions. This brain region is well known for its diverse peptide-expressing cell types, and there are distinct behavioral and physiological consequences of PVH neuropeptide pharmacology as well as cell type–selective optogenetic and chemogenetic perturbations. Molecular and cellular perturbation studies have indicated that the PVH acts as a behavioral state “switchboard” in which each molecularly defined population influences a distinct behavioral state as a labeled-line. However, it is not known if this is consistent with the natural dynamics of PVH neurons during these behaviors.

3We considered three models for the role of molecularly defined cell types to encode behavioral states. The labeled-line model encodes different behavioral states based on the activity pattern of distinct, molecularly defined cell types (Fig. 1A). Labeled-line coding is frequently assumed for information processing in hard-wired neural functions involving essential survival behaviors. At the other end of the spectrum, the full neural ensemble encodes behavioral state irrespective of cell type identity and is primarily a product of neural plasticity (Fig. 1B). In this case, molecular identity is a poor predictor of behavioral state, and functional ensembles are formed using rules that are independent of cell type. An intermediate model is grouped-ensemble coding, in which the members of some molecularly defined cell types show coordinated responses and encode behavioral states based on specific combinations of cell types (Fig. 1C).

4It has been extremely difficult to distinguish these models because there has not been a suitable approach to simultaneously examine the functional dynamics of individual neurons having distinct molecular identities across many behavioral states within a single animal. Several methods have been developed to map gene expression onto neuronal function. However, these methods have limitations of temporal sensitivity, dynamic range, molecular diversity, or restrictions to transparent organisms that have prevented evaluation of this type of problem. Thus, an unbiased and systematic method to relate multigene expression of individual neurons to their activity patterns in vivo on a subsecond time scale is needed to combine molecular and systems neuroscience.

5To address this problem, we developed the CaRMA (calcium and RNA multiplexed activity) imaging platform, a comprehensive, three-part strategy compatible with recording neuron dynamics deep in the mouse brain (Fig. 1D). In the first step, we used single-cell RNA sequencing (scRNA-seq) to classify PVH cells into molecular types (Fig. 1D, a and b). In the second step, we performed volumetric, deep-brain, pan-neuronal, two-photon calcium imaging in awake, behaving animals without regard to PVH neuronal subtypes to measure the response dynamics of hundreds of neurons during multiple behavioral states (Fig. 1D, c). Subsequently, the brain was removed, sectioned, and the neurons in the ex vivo tissue sections were registered to those from the in vivo imaging volume (Fig. 1D, d and e). In the final step, molecular identity of the neurons imaged in vivo was determined post hoc by performing multiple rounds of three-plex RNA-fluorescent in situ hybridization (FISH) on the ex vivo tissue sections, guided by marker genes from scRNA-seq in step 1 (Fig. 1D, f to i). This approach enables monitoring the dynamics of many molecularly defined neurons in parallel, within a deep-brain region, and across multiple behaviors.

6We used scRNA-seq to transcriptionally profile 706 manually picked PVH cells from 10 animals. We applied iterative unsupervised clustering to group cells into 12 molecularly defined cell types within the PVH (fig. S1A). Differentially expressed marker genes for these clusters included canonical PVH neuropeptides: oxytocin (Oxt), vasopressin (Avp), corticotropin-releasing hormone (Crh), thyrotropin-releasing hormone (Trh), somatostatin (Sst), proenkephalin (Penk), and prodynorphin (Pdyn). Within these broad classes, clusters were further divided into subgroups based on the expression of additional marker genes: glutamic acid decarboxylase 2 (Gad2), reelin (Reln), netrin G1 (Ntng1), and neuropeptide Y receptor type 1 (Npy1r) (fig. S1B). All PVH clusters expressed Sim1, contained the excitatory neuron marker Vglut2 (Slc17a6), and lacked the inhibitory neuron marker Vgat (Slc32a1). This targeted, deep-sequencing approach provides a combinatorial set of marker genes for molecularly defined cell type assignment.

7We mapped the spatial and coexpression pattern of these 11 differentially expressed genes plus Vglut2 by performing four rounds of three-plex FISH from two mice in PVH subregions separated by 280 mm: the anterior PVH (aPVH, 1394 cells), middle PVH (mPVH, 1855 cells), and posterior PVH (pPVH, 1263 cells) (Fig. 2A). This was achieved by developing new methods for stripping fluorescently labeled probes, aligning individual cells in three dimensions across multiple rounds of FISH using 4′,6-diamidino-2-phenylindole (DAPI) fluorescence, as well as generalizable methods for three-dimensional (3D) segmentation of cell boundaries (figs. S2 to S4; see the materials and methods).

8We observed marker gene enrichment in different PVH subregions (Fig. 2B): Two genes (Gad2 and Ntng1) were significantly enriched in aPVH, three genes (Vglut2, Crh, and Avp) were significantly enriched in mPVH, and nine genes (Vglut2, Gad2, Npy1r, Crh, Reln, Ntng1, Pdyn, Oxt, and Avp) were significantly enriched in pPVH. Pairwise gene coexpression also differed based on anterior-posterior position (fig. S5). For example, a majority of Crh cells coexpressed Npy1r in pPVH (75%) but considerably less in mPVH (30%) and aPVH (35%).

9We determined the proportions of PVH cells in situ that expressed 0 to 12 of the marker genes (Fig. 2C). Coexpression (more than two genes) was common (81% of cells), and the mode was three genes. Unsupervised hierarchical clustering of gene expression profiles from these 4512 PVH cells resulted in 13 transcriptional clusters (Fig. 2, D to F). The spatial distribution of these molecularly defined cell types was largely intermingled but some well-organized patterns were apparent (Fig. 2G; fig. S6, A and B; and table S1). Expression levels of the marker genes were distributed continuously, with strong rightward skew and high variance (Fig. 2H). Clustering implicitly thresholded gene coexpression relationships, and most FISH clusters (12/13) were dominated by high expression of one gene (Fig. 2E). FISH clusters were correlated with the scRNA-seq dataset (fig. S6, C and D). Likewise, scRNA-seq clusters were represented in the FISH dataset (fig. S6, E to G) and mapped to different regions of the PVH (fig. S6G).

10To record PVH neuronal dynamics, we expressed GCaMP6m in PVH neurons without regard to neuronal subtypes. A thin gradient refractive index (GRIN) lens was implanted into the PVH (Fig. 3A). Volumetric two-photon calcium imaging was acquired from PVH neurons under the GRIN lens (range: 80 to 240 μm) in a head-fixed behaving mouse (Fig. 3B). We performed a set of calcium-imaging experiments involving eating during hunger, drinking water during thirst, hedonic eating (ad libitum fed, palatable food), fear retrieval (movie S1), and hormone-induced hunger (ghrelin) or energy surfeit (leptin) sequentially over 10 days (fig. S7). Subsequently, the brain was removed and sectioned.

11To register the GCaMP-expressing neurons in the ex vivo–sectioned brain to the in vivo image volume, the sections below the GRIN lens were imaged by confocal microscopy. Because the GRIN lens introduced multiple optical distortions, we characterized these aberrations and transformed the in vivo image volume coordinates to the view obtained by confocal imaging in ex vivo tissue sections (fig. S8). For a successful in vivo ↔ ex vivo alignment, it was crucial to correct field curvature introduced by the GRIN lens (fig. S8, C, D, F, and J), as well as the nonlinear relationship between sample position and objective lens position (fig. S8, H and J) and depth-dependent magnification changes (fig. S8, C, D and I). After correcting for these aberrations (Fig. 3C and fig. S9A), we found neurons with distinctive shapes in the two-photon–imaging volume (neurons 1 to 3 in Fig. 3D; neurons 1, 2, 4, and 35 in fig. S9B) that could be recognized in the ex vivo confocal images. We used these neurons as starting points to match the surrounding neurons imaged in vivo and ex vivo with 96% correspondence between two independent observers. Most of the neurons (89.3%, 334/374) imaged in vivo could be found in the brain slices, and 95.5% (319/334) of the aligned neurons could be tracked across the entire sequence of experiments (Fig. 3D and fig. S9B).

12After matching neurons between the in vivo and ex vivo image spaces, we quenched GCaMP6m with acidic buffer and performed four rounds of three-plex RNA-FISH (Fig. 3E and fig. S9C). We developed a non-rigid 3D registration pipeline for the GCaMP-expressing brain slices across multiple rounds of processing (fig. S10). Calcium dynamics of individual neurons during multiple behaviors were extracted from the in vivo image volumes (fig. S11), along with the corresponding molecular profiles for each neuron acquired from ex vivo tissue sections (Fig. 3F and fig. S9D). These computational and image analysis tools are available online (see “Data and materials availability” in the Acknowledgments).

13We recorded calcium dynamics of 319 neurons from three mice in the same portion of the pPVH (fig. S12) during all behaviors. We segmented the behavioral tasks into 11 behavioral states (fig. S7). These behavioral states could be well separated by unsupervised clustering of the recorded PVH neural activity patterns (Fig. 4A).

14First, we investigated whether this PVH ensemble showed evidence for labeled-line coding of behavioral states irrespective of any gene expression information. A labeled-line neuron was defined as specifically activated or inhibited in one state but not in other states. We searched for the maximum number of labeled-line neurons among all k combinations of 11 behavioral states (Fig. 4B). Although selectively tuned neurons were evident with comparison of a small number of behavioral states, the number of neurons showing apparent labeled-line coding was dependent on the number of behavioral states examined (Fig. 4C). Beyond six behavioral statessss, there were no selective neurons for each state. Moreover, the proportions of neurons selective for each behavioral state were statistically greater than chance for not more than four states (Fig. 4B). We also considered the possibility that grouping similar states might make labeled-line coding more apparent. States were grouped into “state classes” either subjectively based on ideas about behavioral similarity (fig. S13A) or systematically based on PVH neuronal response similarity metrics for unsupervised clustering (Fig. 4A). However, grouping further reduced performance of the labeled-line model (Fig. 4D and figs. S13 and S14). Examining a small number of behavioral states might thus lead to the erroneous conclusion that a set of neurons is involved in highly selective labeled-line coding. Although the labeled-line configuration cannot be excluded for some brain functions, labeled-line coding is poorly scalable.

15Next, we used unsupervised clustering of gene expression profiles to group the PVH neurons into 11 molecularly defined clusters (MCs) (Fig. 5, A to C; see the materials and methods) so that we could evaluate the activity patterns in molecularly defined cell types across 11 behavioral states. Although GCaMP expression was not observed in Oxt- or Avp-expressing neurons because of adeno-associated virus tropism (fig. S15), these 11 MCs were consistent with cell types clustered from FISH-only pPVH tissue (fig. S16). Spatially, these molecularly defined cell types were intermingled, primarily in the medial parvicellular and lateral parvicellular subdivisions of the pPVH (Fig. 5, D and E).

16If cell type information is important to encode a behavioral state, then we expected that neurons within a molecularly defined type, grouped solely using gene expression information, should respond similarly in that state. When we looked at MC5-Crh neurons in response to hedonic eating and fear retrieval, we found that neurons within this cluster responded similarly. By contrast, MC8-Penk neurons showed highly heterogeneous responses in these states (Fig. 5F). Thus, molecularly defined neurons do not necessarily respond similarly in all behaviors and may be tuned to only some behavioral states.

17We quantified the similarity of neuronal responses within each cell type cluster using linear purity (fig. S17; see the materials and methods). During fear cue presentation, MC5-Crh had high purity (approaching 1), whereas MC8-Penk had low purity (close to zero; Fig. 5F). Cell types with purity >0.5 can be considered to respond similarly (fig. S17). Next, we defined a metric called “consistent-response”, which weights cell type activity by the response purity (Fig. 5F and fig. S17). A high consistent-response molecularly defined cell type is composed of neurons that respond strongly to a stimulus with similar temporal dynamics. Low consistent-response cell types have either low purity (a diversity of responses) or low response magnitude in a behavioral state (fig. S17).

18We computed the temporal maximum of consistent-responses and purities of each cell type in 11 behavioral states (Fig. 5H and fig. S18, B and C). In the consistent-response map, some cell types responded differently to individual behavioral states, whereas others were similarly tuned across multiple states. MC8-Penk neurons, which had relatively low response and low purity during fear, showed higher purity activation after injection of the hormone ghrelin. This demonstrates that molecularly defined neurons are tuned to respond similarly in some behaviors but not others. The purity of the molecularly defined cell types resulting from unsupervised clustering of multiple genes produces significantly higher response purities than when cells are grouped by expression of a single gene (fig. S19, A to C). The high purity activation response to food consumption of MC5-Crh neurons that highly express Crh and coexpress Vglut2 and Npy1r in the pPVH is opposite to the inhibitory response reported by low-resolution, multineuron fiber photometry measurements of all CRH neurons in the PVH, although both methods give the same response to fear-inducing stimuli. However, low-expressing Crh+/Vglut2+/Npy1r− neurons in the pPVH were inhibited in response to food ingestion (fig. S20), and this cell type is more prevalent in the mPVH (fig. S20, F and G), which is the major target region in the reported photometry experiments. Thus, CaRMA imaging is a powerful platform for systematically discovering the relationship between cell types defined by multiple gene markers and their response consistency in different behavioral states.

19The consistent-response map shows the relationships between molecularly defined cell type activity patterns for behavioral state coding (Fig. 5H). From these data, it is apparent that behavioral state is not represented by a single molecularly defined cell type (fig. S19D). Instead, the consistent-response pattern indicates combinatorial coding. What is different from previous views of hypothalamic circuits is that behavioral state is encoded by the different combinations, magnitudes, and temporal dynamics of cell type activity patterns. For example, MC5-Crh and MC6-Pdyn showed similar consistent-responses during hunger eating, drinking, and hedonic eating (absolute differences ≤0.18; Fig. 5H), but their temporal response profiles were significantly different (Fig. 5G). MC5-Crh neurons ramped to maximum magnitude, whereas MC6-Pdyn neurons responded more quickly to the ingested stimuli (Fig. 5G). Similarly, the return to baseline after withdrawing food, water, or the offset of the fear retrieval cue was slower for MC5-Crh (fig. S18A). Although these cell types have similar consistent-response measures, their different temporal dynamics indicate distinct functions for encoding behavioral states, as is also the case for other cell types (fig. S18A).

20CaRMA imaging enabled us to investigate whether behavioral states can be quantitatively decoded from the PVH ensemble dynamics of molecularly defined cell types. Because cell type encoding of behavioral states has been traditionally examined either one cell type at a time or was used without regard to cell type, we first examined the decoding performance of individual molecularly defined cell types to distinguish 11 behavioral states.

21We used multinomial logistic regression classification to determine the overall performance of predicting behavioral states from the temporal dynamics of single neurons without regard to cell type (All) as well as from singlssse neurons of known cell types (Fig. 6A; see the materials and methods). Single neurons from several cell types showed decoding performance superior to a classifier trained using all neurons without regard to cell type (Fig. 6B and fig. S21A). The decoding accuracy from the dynamics of MC5-Crh neurons was the highest (Fig. 6B), and the confusion matrix for behavioral state decoded by MC5-Crh neurons showed perfect decoding of the ghrelin-induced hunger state and high decoding performance for fear retrieval (Fig. 6C). MC5-Crh neurons showed significantly different response amplitude and sign across several behavioral states (Fig. 6D), which explained the high decoding performance. However, our analysis also showed that individual cell types in the PVH are insufficient to decode many behavioral states.

22CaRMA imaging enables the activity of many cell types to be imaged simultaneously in the same animal. This permitted us to test models of behavioral state coding using an ensemble of molecularly defined cell types (Fig. 1C) or, alternatively, an ensemble of neurons in which molecularly defined cell type information was ignored (Fig. 1B). We trained a classifier on an ensemble of molecularly defined cell types model by randomly choosing one cell’s activity trace from each of 10 molecularly defined cell types and concatenating them into 10-cell ensembles for each behavioral state (repeated 100 times with replacement), followed by multinomial logistic regression for 11 behavioral states using 10-fold cross-validation (accuracy: 93.13 ± 0.02%, Fig. 6, E and F; see the materials and methods). With molecularly defined cell type information, 11 different behavioral states could be accurately decoded with the combinatorial responses of only 10 PVH neurons, one from each of the 10 cell types. This was significantly higher decoding accuracy (P = 0; Fig. 6G) than a neural ensemble that ignores molecularly defined cell type information by randomly assigning neurons into dummy cell types (fig. S21B).

23We also examined the effect of adjusting the hierarchical molecularly defined cell type clustering threshold on behavioral state decoding. Purity and decoding accuracy generally increased with the number of molecular clusters (figs. S22 and S23). Although grouped-ensemble coding was supported with a range of clustering thresholds, this approach involves “peeking” at the functional responses to supervise molecular clustering. Because our objective was to evaluate the predictive value of molecular information on functional and behavioral responses, we proceeded using the unbiased determination of 11 MCs in the pPVH ensemble (to avoid circular logic). Taken together, our findings support a grouped-ensemble coding strategy used by molecularly defined PVH neurons (Fig. 1C).

24To reveal the relative contributions of different molecularly defined cell types to grouped-ensemble coding (Fig. 1C) for different behavioral states in the PVH, we used the multinomial regression coefficients and the consistent-response measurements to generate cell type response-decoding diagrams for each behavioral state (Fig. 6, H to K, and fig. S24). In these diagrams, cell type consistent-responses are weighted by each cell type’s contribution to decoding accuracy and quantitatively summarized across multiple behavioral states. MC5-Crh and MC6-Pdyn neurons had large consistent-responses in several behavioral states. These neurons were activated by negative and positive stimuli, indicating an unexpected general role for encoding salience of diverse stimuli in a variety of different states. Other molecularly defined cell types, such as MC8-Penk and MC10-Trh neurons, contributed substantial weights for distinguishing behavioral states (Fig. 6, H to K, and fig. S24C). For example, food withdrawal (Eat-Aft) in hunger elicited a transient inhibitory response in MC10-Trh after spout retraction that transitions to a prolonged activated response, which was not observed after removing water in thirst (Drink-Aft) (fig. S25A). This contributes to behavioral state coding (fig. S25B) and is also consistent with past work indicating that TRH neurons project to the arcuate nucleus and activate hunger-associated Agouti-related peptide (AGRP) neurons. By contrast, most MC8-Penk neurons showed prolonged inhibition after food withdrawal in hunger but only transient inhibition after drinking in thirst (fig. S25, C and D). More generally, these response-decoding diagrams show a functional organization in the pPVH, with major branches associated with the sign and amplitude of broadly tuned MC5-Crh and MC6-Pdyn neurons and individual twigs associated with a more selective combinatorial code of cell types that distinguish behavioral states.

25Grouped-ensemble coding relies on similar responses from the neurons that make up a cell type, thereby simplifying the neural ensemble to component cell types instead of the much larger total number of cells. Consistent-response maps also indicate higher-level organization based on similar responses between groups of cell types. For example, MC5-Crh and MC6-Pdyn show similar response tuning to multiple behavioral states. By contrast, MC8-Penk typically responds oppositely to these cell types (Figs. 5H and 6, H to K, and fig. S24C). Thus, different molecularly defined cell types may have coordinated activity patterns that comprise hierarchical functional groupings in the neural ensemble. This raises the second major question: Do differentially expressed genes map to functional response types in the neural ensemble?

26We aimed to identify optimal combinations of genes that best predict the functional response classes within the neuronal ensemble. Moreover, neurons have complex dynamics during behaviors, so we assessed the temporal window over which molecular markers offer useful predictive information about functional responses.

27For each behavioral state, we grouped all neurons on the basis of their neuronal response dynamics into functional clusters (FCs) by unsupervised hierarchical clustering (Fig. 7, A and B, and fig. S26). Two FCs emerged for each behavioral state associated with primarily activated (Act-FC) or inhibited (Inh-FC) neurons. We first examined how functional clustering segregates neuron gene expression profiles.

28A subset of gene expression profiles specifically segregated with functional response clusters (Fig. 7, A to C, and fig. S26). We performed differential gene expression analyses between the two FCs in each behavioral state (fig. S27A), which was measured by either the P value of the expression-level difference or by the ratio of expression level between FCs (Fig. 7D and fig. S27B). Npy1r was significantly enriched across 11 behavioral states and was ranked as the first or second enriched gene by P value. Crh and Pdyn were significantly enriched in eight and nine behavioral states, respectively. Moreover, if Npy1r, Pdyn, and Crh were enriched for a behavioral state, then they were found in the same FC because of partial coexpression in the same cells. For example, Npy1r, Pdyn, and Crh were enriched in the Act-FC for Eat-Hed (hedonic) and Fear states. By contrast, the Inh-FC was enriched with neurons expressing Penk. We also examined combinations of molecularly defined cell types within FCs in each behavioral state. For most behavioral states (eight of 11 states), activated or inhibited FCs also showed significant enrichment of cell types (fig. S28). The cell types with the highest FC selectivity were MC5-Crh and MC6-Pdyn, which also showed high consistent-responses in most behavioral states (Fig. 5H). These gene and cell type enrichment analyses demonstrated a strong association between neuron function during behavior and the molecular profile of individual neurons. However, this analysis lacks useful metrics for ranking the predictive power of gene expression for functional response class, and it does not define the optimal set of marker genes for predicting responses across the neural ensemble.

29We therefore determined the best gene combinations for predicting neuron functional response classes in each behavioral state (Fig. 8A). We used supervised machine learning to predict the FCs of individual PVH neurons solely by their gene expression profiles with a logistic regression classifier (Fig. 8B; see the materials and methods). Three normalized features represented the expression level of each gene (figs. S29 and S30A). This approach does not binarize gene expression by thresholding. Instead, it explicitly incorporates the broad continuous distribution of gene expression levels that we found (Fig. 2H) as the independent variables in the predictive model. Sequential forward feature selection (SFFS) determines the most predictive individual gene for FC classification and sequentially adds the next most predictive gene until predictive accuracy no longer increases. This approach identified optimal gene sets that best predicted the functional response class while minimizing the contribution of redundant gene coexpression information (Fig. 8C and fig. S30, B to D). We compared the predictive power of the molecular profiles for functional response class with the 95th percentile of predictive accuracies after shuffling the relationship between gene expression of individual neurons and their FCs (Fig. 8, C and D, and fig. S30B). The two behavioral states during which neuronal ensemble response could be best predicted by gene expression profiles were during ghrelin-induced hunger [accuracy: 73%, area under the receiver-operating characteristic (auROC): 0.784, P = 0) and fear retrieval (accuracy: 69.9%, auROC: 0.745, P = 0) states. The lowest predictive accuracy was for saline injection (accuracy: 62.1%, auROC: 0.570, P = 0.0005) (Fig. 8D and fig. S30D).

30Combinations of genes always had higher predictive power for functional response class than single genes (fig. S30B). In addition, for all behavioral states, the optimal accuracies from the first round of SFFS (SFFS1) with just two to six genes were higher than the accuracies predicted with all genes (Fig. 8C). This showed that a subset of genes provided higher prediction performance than the whole set of genes. Next, we quantitatively assessed and ranked the predictive accuracy of single genes for FCs across 11 behavioral states. Pdyn, Npy1r, or Crh showed the highest predictive power in the first round for all behavioral states. Npy1r was always in the optimal gene set from SFFS1 in all behavioral states. Moreover, Npy1r expression strongly predicted fear-activated neurons. We determined the roles of individual genes within the optimal gene set to improve prediction performance by measuring the average coefficients of those genes from logistic regression (Fig. 8C). For ghrelin-induced hunger, SFFS1 showed an additive relationship for Npy1r and Pdyn prediction of the inhibited functional response class (Fig. 8C). Both genes have been previously associated with appetite control by the PVH. For food consumption during hunger (Eat-Hng), Pdyn and Npy1r predicted the activated response type, whereas the inhibited response class was predicted by Trh and Penk (Fig. 8C). For fear retrieval, Npy1r and Crh best predicted activated neurons and Penk predicted fear-inhibited neurons (Fig. 8C).

31Multiple rounds of SFFS (mSFFS; see the materials and methods) revealed genes that could compensate for the removal of the most predictive gene from the previous rounds (fig. S30B). This indicated that the most predictive genes could be largely compensated by combinations of additional genes in subsequent rounds. For example, Pdyn was most predictive during Eat-Hng, but this predictive accuracy was largely replaced by Crh in round 2 and Npy1r in round 3 of mSFFS, and we found that Npy1r neurons coexpressed Pdyn and Crh in a hierarchical relationship in pPVH (fig. S30, E and F). Consistent with this, Npy1r showed statistically significant predictive power for every behavioral state, but if Npy1r was removed from the analysis, its predictive power could be largely replaced by other genes (fig. S30B). Thus, the hierarchical gene expression relationship of these three genes in the pPVH was associated with similar neuron functional responses. Neurons with low or no Npy1r expression typically responded oppositely. We also found similar results for the predictive power of marker genes for neuron FCs concatenated across all 11 behavioral states (fig. S31).

32Using the quantitative gene expression level to predict functional responses also offers a systematic method for defining the gene expression threshold that targets a functional response class with the best precision for each behavioral state. We identified optimal cut points for the gene expression threshold using Youden’s J statistic with the ROC curve (fig. S32). Prediction precision with a single gene was up to 40% higher, and response purity was significantly improved when the optimal expression threshold was used instead of a lower threshold based on expression over background (Fig. 8, E to G). Thus, neuronal functional identity is not a binary quality of gene expression.

33Next, we investigated how the predictive power of gene expression profiles was related to neuron temporal dynamics. For example, does a gene have similar predictive value throughout a state? For each timestamp (0.4-s interval), we modeled the relationship of gene expression level to neuronal response (Fig. 9, A to C, and fig. S33A). Our results complemented those obtained from prediction of FCs, including (i) higher predictive power of combinatorial expression profile (Fig. 9A), (ii) improvement in prediction performance by mSFFS (Fig. 9A), (iii) redundant information in molecular profiles for predicting temporal response (Fig. 9B), and (iv) that Npy1r predictive power was high and statistically significant across all behavioral states (Fig. 9A). However, predictive power depended on cell type–specific neuronal dynamics. During eating and drinking, the predictive power of marker genes increased as consumption progressed. For example, Pdyn neuron activity rose more quickly than Pdyn predictive accuracy for Eat-Hed and Eat-Hng. In this case, the transient activation at onset of food presentation was not specific within Pdyn-expressing neurons but was also widespread in Pdyn-negative neurons (fig. S33B). This shows that cell type–specific responses develop progressively in some behaviors. By contrast, fear retrieval showed an immediate cell type–specific response that was best predicted by Npy1r and gradually became less well predicted by gene expression (fig. S33B).

34In five of 11 behavioral states (Eat-Hng, Drink-Thirst, Eat-Hed, Fear, and Ghrelin), a single gene had the greatest predictive accuracy throughout a behavioral state. However, for six other states, the most predictive gene for the ensemble functional response switched with time. Moreover, using the optimal predictions of neural activity, we ranked neurons by the fraction of deviance explained with our linear regression models along the entire time series across behavioral states (Fig. 9, C and D, fig. S33A). A majority of neurons (65%) showed statistically significant predictive power. Neurons within the top quartile of predictive powers showed significantly higher expression levels of Vglut2, Gad2, Npy1r, Crh, Reln, Pdyn, and Trh (fig. S33D). Next, we remapped the corresponding gene expression profiles onto molecularly defined cell types (Fig. 9, E and F). MC4-Npy1r, MC5-Crh (exclusively), and MC6-Pdyn neurons were enriched in the most highly predictive quartile (fig. S33E), consistent with the high predictive accuracy associated with Npy1r coexpression.

35Neural coding by an ensemble-of-cell-types indicates that many molecularly defined neurons also form functional groupings that work in concert to encode behavioral state. It is unlikely that coding in the pPVH follows the labeled-line model because we could not find such neurons in the ensemble irrespective of molecular identity. Our results validate the use of molecularly defined cell types to reduce the dimensionality of the PVH neuronal ensemble to understand behavioral state coding. Ensemble-of-cell-types coding of behavioral states provides superior efficiency, scalability, and flexibility than labeled-line coding and offers superior robustness (due to redundancy) than cell type–independent, full-ensemble coding. Operationally, grouping by molecularly defined cell types is justified by high consistent-responses and is analogous to the pooling operation in artificial neural networks, which is used to extract invariant information. This may also be a key function when molecularly defined circuits are used for coding behavioral states, which should be invariant to irrelevant external conditions. Taken together, these cell types and the patterns of activity reported here provide a parts list for systematically examining the neural coding of behavioral state by the PVH. CaRMA imaging offers a quantitative method to evaluate the extent to which other brain regions rely on molecularly defined cell type ensemble coding.

36Using CaRMA imaging, we demonstrate an important role for gene expression information to predict both neuron functional response class as well as temporal dynamics within the PVH neural ensemble. High predictive power of Npy1r for neuron activity was observed with all analysis methods. However, combinations of other cell type marker genes also showed significant prediction for functional responses. The highest consistent-response cell types, MC5-Crh and MC6-Pdyn, coexpressed Npy1r, and cell types that responded oppositely typically lacked Npy1r. scRNA-seq analysis revealed hierarchical relationships of gene expression within neuron populations, where the smallest subpopulations are designated as “types” and the largest subdivisions as “classes” that encompass many types. We found evidence that this molecular hierarchy extends to functional responses in the neuronal ensemble, in which similarly tuned neuron types form classes that can be well described by a neuromodulator receptor gene that is expressed across multiple molecularly defined cell types (fig. S30, E and F). This is interesting in the context of the neuromodulatory function of Npy1r and its endogenous agonist NPY. Because Npy1r is expressed across multiple cell types and NPY signals through volume transmission, this receptor can cooperatively regulate diverse cell types (Fig. 10A). This may enhance the generality of neural output in the presence of noisy neural inputs, which is important to facilitate behavioral state coding but is also a computational role akin to regularization in machine learning. Nevertheless, despite the response similarity between neurons that express Npy1r, there are reliable temporal differences that map onto cells with specific gene expression profiles. These temporal responses are important for distinguishing behavioral states and functional response types in the neuronal ensemble. In the analogy of an orchestral ensemble, Npy1r is a conductor selectively coordinating the responses of multiple cell types (Fig. 10B) that are like instruments that play together but maintain distinctive contributions to a score.

37Marker-gene expression in the PVH can be dynamic, but to the extent that this affects our results, it may degrade the predictive accuracy of gene expression level to a functional or behavioral response. Additional genetic markers may improve these predictions by offering increased coverage of all neurons in the ensemble. However, predictive accuracies solely from gene expression information are not expected to be perfect because functional responsiveness is also influenced by plasticity and connectivity mechanisms that are independent of cell type. Moreover, our findings also highlight potential limitations of site-specific recombinase (e.g., Cre) mouse lines that use low gene expression thresholds to specify cell type identity.

38Although we focused here on decoding behavioral states and predicting functional responses solely from molecular information about the underlying cells, we found that if functional information is used to supervise clustering, then molecularly defined cell types can be defined with higher purity. This may be useful in the future for better defining combinations of molecular markers and their expression levels that are optimal for targeting specific functional perturbations.

39CaRMA imaging provides a solution to the general problem of assessing functional roles for the large number of new molecularly defined cell types that are being discovered routinely using single-cell transcriptomics. Because CaRMA imaging does not rely on genetic modifications, only viral targeting, it is applicable to any organism suitable for in vivo functional imaging. By producing unbiased and systematic datasets that integrate genes, neuron dynamics, and behavior, CaRMA imaging offers a path to discoveries about the organization of brain function at both the molecular level and the systems level.