Introduction

Individuality is a fundamental aspect of behavior that is observed even among genetically-identical animals reared in similar environments. We are specifically interested in individuality that is evident as idiosyncratic differences in behavior that persist for much of an animal’s lifespan. Such variability is observed across species including round worms (Stern et al., 2017), aphids (Schuett et al., 2011), fish (Laskowski et al., 2022), mice (Freund et al., 2013), and people (Johnson et al., 2010). Small, genetically tractable model species, such as Drosophila, are particularly promising for discovering the genetic and neural circuit basis of individual behavior variation. Flies exhibit individuality in many behaviors (Werkhoven et al., 2021), and the mechanistic origins of this variation has been studied for phototactic preference (Kain et al., 2012), temperature preference (Kain et al., 2015), locomotor handedness (Ayroles et al., 2015; Buchanan et al., 2015; de Bivort et al., 2022), object-fixated walking (Linneweber et al., 2020), and odor preference (Honegger et al., 2019). Generally, the neural substrates of individuality are poorly understood, though in a small number of instances nanoscale circuit correlates of individual behavioral biases have been identified (Lillvis et al., 2022; Linneweber et al., 2020; Skutt-Kakaria et al., 2019). We hypothesize that as sensory cues are encoded and transformed to produce motor outputs, their representation in the nervous system becomes increasingly idiosyncratic and predictive of individual behavioral responses. We seek to identify “loci of individuality” – sites at which this idiosyncrasy emerges.

Olfaction in the fruit fly Drosophila melanogaster is an amenable sensory system for identifying loci of individuality, as 1) individual odor preferences can be recorded readily, 2) neural representations of odors can be measured via calcium imaging, 3) the circuit elements of the pathway are well-established, and 4) a deep genetic toolkit enables mechanism-probing experiments. The neuroanatomy of the olfactory system, from the antenna through its first central-brain processing neuropil, the antennal lobe (AL), is broadly stereotyped across individuals (Couto et al., 2005; Grabe et al., 2015; Wilson et al., 2004). The AL features ∼50 anatomically identifiable microcircuits called glomeruli (Figure 1A). Each glomerulus represents an odor-coding channel and receives axon inputs from olfactory receptor neurons (ORNs) expressing the same olfactory receptor gene (de Bruyne et al., 2001). Uniglomerular projection neurons (PNs) carry odor information from each glomerulus deeper into the brain (Jeanne and Wilson, 2015). AL-intrinsic local neurons (LNs) project among glomeruli (Chou et al., 2010) and modulate odor representations (Wilson and Laurent, 2005). Glomerular organization is a key stereotype of the AL; using glomeruli as landmarks, one can identify comparable ORN axons and PNs across individuals.

Figure 1.

Individual flies differ in their PN calcium responses to identical odor stimuli, as well as their odor-vs-odor preference choices (Honegger et al., 2019). Several possible determinants of individual odor preference can already be hypothesized for the fly olfactory circuit (Rihani and Sachse, 2022). The extent of preference variability depends on dopamine and serotonergic modulation (Honegger et al., 2019). Neuromodulation clearly plays a role in the regulation of behavioral individuality (Maloney, 2021), but its effects vary by modulator and behavior (de Bivort et al., 2022; Kain et al., 2012). With respect to wiring variation, the number of ORNs and PNs innervating a given glomerulus varies within hemispheres (Tobin et al., 2017) and across individuals (Grabe et al., 2016; Schlegel et al., 2020), as does the glomerulus-innervation pattern of individual LNs (Chou et al., 2010). Subpopulations of LNs and PNs express variable serotonin receptors (Sizemore and Dacks, 2016), so the effects of neuromodulation and wiring may interact to influence individuality. Little is known about possible molecular or nanoscale correlates of individual behavioral bias. Thus, individual odor preference could have its origins in many potential mechanisms, ranging from circuit wiring to modulation to neuronal intrinsic properties.

Outside the olfactory system, there are a few examples in which microscale circuit variation predicts individual behavioral preference. Wiring asymmetry in an individual fly’s dorsal cluster neurons is predictive of the straightness of its object-oriented walking behavior (Linneweber et al., 2020), and left-right asymmetry in the density of presynaptic sites of protocerebral bridge to lateral accessory lobe-projecting neurons predicts an individual fly’s idiosyncratic turning bias (Skutt-Kakaria et al., 2019). The number of synaptic connections from the pC2l to pIP10 neurons correlates with male song rate during courtship (Lillvis et al., 2022), and the presence of ectopic branches in neurons of the T2 hemilineage predicts delayed spontaneous flight initiation (Mellert et al., 2016).

In this work, we sought to identify loci of individuality by measuring odor preferences and neural responses to odors in the same individuals and determining the extent to which the latter predicted the former. We found that idiosyncratic calcium responses in PNs were correlated with individual preferences in a choice between two aversive odorants. Examining a molecular component presynaptic to PNs, we found that the density of the scaffolding protein Bruchpilot also predicts odor preference. To unify these results and connect wiring variation to circuit outputs and behavior, we simulated developmental variation in a 3,062-neuron spiking model of the antennal lobe. Simulated stochasticity in the ORN-PN synapse recapitulated our empirical findings. Thus, we identified the ORN-PN synapse as a likely locus of individuality in fly odor preference, demonstrating that behaviorally-relevant variation in neural circuits can be found in the sensory periphery at the nanoscale.

Results

Individual flies encode odors idiosyncratically

Focusing on behavioral variation within a genotype, we used isogenic animals expressing the fluorescent calcium reporter GCamp6m (Chen et al., 2013) in either of the two most peripheral neural subpopulations of the Drosophila olfactory circuit, ORNs or PNs (Figure 1E). We performed head-fixed 2-photon calcium imaging after measuring odor preference in an untethered assay (Honegger et al., 2019) (Figure 1B-D, Figure 1 – figure supplement 1A). Individual odor preferences are stable over timescales longer than this experiment (Figure 1 – figure supplement 1B-E).

We measured volumetric calcium responses in the antennal lobe (AL), where ORNs synapse onto PNs in ∼50 discrete microcircuits called glomeruli (Figure 1A) (Couto et al., 2005; Grabe et al., 2015). Flies were stimulated with a panel of 12 odors plus air (Figure 1D, Figure 1 – figure supplement 2) and k-means clustering was used to automatically segment the voxels of 5 glomeruli from the resulting 4-D calcium image stacks (Figure 1E, Figure 1 – figure supplement 5, Materials and Methods) (Couto et al., 2005). Both ORN and PN odor responses were roughly stereotyped across individuals (Figure 1G,H), but also idiosyncratic (Honegger et al., 2019).

Figure 2.

Responses in PNs appeared to be more idiosyncratic than ORNs (Figure 1J); a logistic linear classifier decoding fly identity from glomerular responses was more accurate when trained on PN than ORN responses (Figure 1 – figure supplement 6A). While the responses of single ORNs are known to vary more than those of single PNs (Wilson, 2013), our recordings capture the total response of all ORNs or PNs in a glomerulus. This might explain our observation that ORNs exhibited less idiosyncrasy than PNs. PN responses were more variable within flies, as measured across the left and right hemisphere ALs, compared to ORN responses (Figure 1 – figure supplement 6C; p < 2x10^-5, Mann-Whitney U test), suggesting that odor representations become more divergent farther from the sensory periphery.

PN, but not ORN, responses predict odor-vs-odor preference

Next we analyzed the relationship of idiosyncratic coding to odor preference, by asking in which neurons (if any) did calcium responses predict individual preferences of flies choosing between two aversive monomolecular odors: 3-octanol (OCT) and 4-methylcyclohexanol (MCH).

Because we could potentially predict preference (a single value) using numerous glomerular-odor predictors, and had a limited number of observations (dozens), we used dimensionality reduction to hold down the number of comparisons we made. We computed the principal components (PCs) of the glomerulus-odor responses (in either ORNs or PNs) across individuals (Figure 1G-I; Figure 1 – figure supplement 3, Figure 1 – figure supplement 8) and fit linear models to predict the behavior of individual flies from their values on the odor response PCs. No PCs of ORN neural activity could linearly predict OCT-MCH preference beyond the level of shuffled controls (n = 35 flies) (Figure 1K,L). The best ORN PC model only predicted odor-vs-odor behavior with a nominal R² of 0.031. In contrast, PC 2 of PN activity was a statistically significant predictor of odor preference, accounting for 15% of preference variance in a training set of 47 flies (p = 0.0063; Figure 2 – figure supplement 1C) and 31% of preference variance on test data of flies (p = 0.0069; Figure 2 – figure supplement 1D). These p-values remaining significant at 𝛼 = 0.05 following a Bonferroni correction for 5 comparisons. Combined train/test statistics (r = 0.45; p = 0.0001) are presented in Figure 1K,M. Thus, idiosyncratic PN calcium predicts odor vs. odor preference.

Figure 3.

We conducted a follow-up analysis to contextualize the finding of calcium PCs predicting odor preference with r = ∼0.4. This value is lower than 1.0 due to at least two factors: 1) any non-linearity in the relationship between calcium responses and behavior, and 2) sampling error in, and temporal instability of, behavior and calcium responses over the duration of the experiment. A lower bound on the latter can be estimated from the repeatability of behavioral measures over time (Figure 1 – figure supplement 1B-E). We performed a statistical analysis to roughly estimate model performance if there were no sampling error or drift in the measurement of behavior and calcium responses (Figure 1 – figure supplement 9; Materials and Methods). This analysis suggests that the measured correlation between calcium and behavior (⍴_signal) would be 0.68 in the absence of sampling error and temporal instability, but the uncertainty in this estimate is high (90% CI: 0.24-0.95).

We additionally assessed the extent to which idiosyncratic calcium responses in ORNs or PNs could predict preference between air and a single aversive odor (OCT). We found a suggestive correlate: PC 1 of ORN calcium responses explained 23% of preference variance (n = 30 flies, p = 0.0099, Figure 1 – figure supplement 10B), but this association was dominated by a single outlier (R² of 0.078, p = 0.14 with the outlier removed).

We next sought a biological understanding of the models associating calcium responses with odor preference. The loadings of the ORN and PN PCs indicate that variation across individuals was correlated at the level of glomeruli much more strongly than odorant (Figure 1H; Figure 1 – figure supplements 3, 8). This suggests that stochastic variation in the olfactory circuit results in individual-level fluctuations in the responses of glomeruli-specific rather than odor-specific responses. In the odor-vs-odor preference model, the loadings of PC2 of PN calcium responses contrast the responses of the DM2 and DC2 glomeruli with opposing weights (Figure 2A), suggesting that the activation of DM2 relative to DC2 predicts the likelihood of a fly preferring OCT to MCH. Indeed, a linear model constructed from the average DM2 minus average DC2 PN response (Figure 2B) showed a statistically significant correlation with preference for OCT versus MCH (R² = 0.12; p = 0.0035; Figure 2C). The model slope coefficient was negative (Table 1), indicating that greater activation of DM2 vs DC2 correlates with preference for MCH. With respect to odor-vs-odor behavior, we conclude that the relative responses of DM2 vs DC2 in PNs compactly predict an individual’s preference.

Table 1:

Odor experience has been shown to modulate subsequent AL responses (Golovin and Broadie, 2016; Iyengar et al., 2010; Sachse et al., 2007). This raises the possibility that our models were actually predicting individual flies’ past odor experiences (i.e., the specific pattern of odor stimulation flies received in the behavioral assay) rather than their preferences. To address this, we imposed the specific odor experiences of previously tracked flies (in the odor-vs-odor assay) on naive “yoked” control flies (Figure 2D) and measured PN odor responses of the yoked flies. Applying the PN PC 2 model to the yoked calcium responses did not predict flies’ odor experience (R² = 0.019, p = 0.49; Figure 2E). This is consistent with PN calcium responses predicting odor preference rather than odor experience.

(Mazor and Laurent, 2005) found that PN response transients, rather than fixed points, contain more odor identity information. We therefore asked at which times during odor presentation an individual’s neural responses could best predict odor preference. Applying our calcium-to-behavior models (PN PC2-odor-vs-odor, as well as ORN PC1-odor-vs-air, PN PC1-odor-vs-air) to the time-varying calcium signals, we found that in all cases, behavior prediction rose during odor delivery (Figure 2 – figure supplement 2). In ORNs, the predictive accuracy remained high after odor offset, whereas in PNs it declined. The times during which calcium responses predicted individual behavior generally aligned to the times during which a linear classifier could decode odor identity from neuronal responses (Figure 2 – figure supplement 2D), suggesting that idiosyncrasies in odor encoding predict individual preferences.

Variation in a presynaptic scaffolding protein predicts odor preference

We next investigated how structural variation in the nervous system might relate to idiosyncratic behavior. Because PN, but not ORN, calcium responses predicted odor-vs-odor preference, we hypothesized that a circuit element between ORNs to PNs could confer onto PNs behaviorally-relevant physiological idiosyncrasies absent in ORNs. We therefore imaged presynaptic T-bar density in ORNs using transgenic mStrawberry-tagged Brp-Short, immunohistochemistry and confocal microscopy (Mosca and Luo, 2014) after measuring individual preference for OCT versus MCH (Figure 3A). Brp-Short density was quantified as total fluorescence intensity / glomerulus volume for 4 of the 5 focus glomeruli (Figure 3B, Figure 3 – figure supplement 1A-F; DL5 was not readily segmentable in our confocal samples). We chose this metric as we found it could be used to predict individual behavioral biases in a previous study (Skutt-Kakaria et al., 2019). This measure was consistent across hemispheres (Figure 3 – figure supplement 1C), while also showing variation among individuals, like calcium responses.

To relate presynaptic structural variation and behavior, we used the same analytical approach as we had for calcium responses. PCs 1 and 2 of Brp-Short density had notable similarities to those of the calcium responses: PC 1 was positive across glomeruli and PC 2 exhibited a sign contrast between DC2 loadings and all other glomerulus loadings (Figure 3 – figure supplement 1G). As in the PN calcium response models, PC 2 of Brp-Short density was the best predictor of odor-vs-odor preferences in training data (Figure 3D-E, Figure 3 – figure supplement 1I, R² = 0.22, n = 22 flies, p = 0.028) and for test data (Figure 3 – figure supplement 1J, R² = 0.078, n = 31 flies, p = 0.13; statistics from combined train and test data: R² = 0.088, n = 53 flies, p = 0.031; Figure 3F; median ⍴_signal 0.39, 90% CI 0.03-0.86). To better understand the microstructural basis of our Brp-Short density metric, we performed paired behavior and expansion microscopy (Asano et al., 2018; Gao et al., 2019) in flies expressing Brp-Short specifically in DC2-projecting ORNs (Supplementary Video 3). Expansion yielded a ∼4-fold increase in linear resolution, allowing imaging of individual Brp-Short puncta (Figure 3 – figure supplement 1K). While the sample size (n = 8) of this imaging pipeline was insufficient for a formal statistical analysis, the trend between Brp-Short density in DC2 (measured as individual puncta / glomerular volume) and odor-vs-odor preference was more consistent with a positive correlation than other metrics, such as median puncta volume (Figure 3 – figure supplement L,M).

The best presynaptic density models are less predictive of behavior than the best calcium response models (R² = 0.088 vs R² = 0.22; ⍴_signal = 0.39 and 0.68, respectively; Figure 2 – figure supplement 1C,D vs Figure 3 – figure supplement 1I,J), suggesting that presynaptic density variation is not the full explanation of calcium response variability. Nevertheless, differences in presynaptic inputs to PNs may contribute to variation in the calcium dynamics of those neurons, in turn giving rise to individual preferences for OCT versus MCH.

Developmental stochasticity in a simulated AL recapitulates empirical PN response variation

Finally, we sought an integrative understanding of how synaptic variation plays out across the olfactory circuit to produce behaviorally-relevant physiological variation. We developed a leaky-integrate-and-fire model of the entire AL, comprising 3,062 spiking neurons and synaptic connectivity taken directly from the Drosophila hemibrain connectome (Scheffer et al., 2020).

After tuning the model to perform canonical AL computations, we introduced different kinds of stochastic variations to the circuit and determined which (if any) would produce the patterns of idiosyncratic PN response variation observed in our calcium imaging experiments (Figure 4A). This approach assesses potential mechanisms linking developmental variation in synapses to physiological variation that apparently drives behavioral individuality.

Figure 4.

The biophysical properties of neurons in our model (Figure 4B, Table 2) were determined by published electrophysiological studies (See Voltage model in Materials and Methods) and were similar to those used in previous fly models (Kakaria and de Bivort, 2017; Pisokas et al., 2020). The polarity of neurons was determined largely by their cell type (ORNs are excitatory, PNs predominantly excitatory, and LNs predominantly inhibitory – explained further in Materials and Methods). The strength of synaptic connections between any pair of AL neurons was given by the hemibrain connectome (Scheffer et al., 2020) (Figure 4C). Odor inputs were simulated by injecting current into ORNs to produce spikes in those neurons at rates that match published ORN-odor recordings (Münch and Galizia, 2016), and the output of the system was recorded as the firing rates of PNs during odor stimulation (Figure 4D). At this point, there remained only four free parameters in our model, the relative sensitivity (postsynaptic current per upstream action potential) of each AL cell type (ORNs, PNs, excitatory LNs and inhibitory LNs). We explored this parameter space manually, and identified a configuration in which AL simulation (Figure 4 – figure supplement 1) recapitulated four canonical properties seen experimentally (Figure 4 – figure supplement 2): 1) typical firing rates at baseline and during odor stimulation (Bhandawat et al., 2007; Dubin and Harris, 1997; Jeanne and Wilson, 2015; Seki et al., 2010), 2) a more uniform distribution of PN firing rates compared to ORN rates (Bhandawat et al., 2007), 3) greater separation of PN odor representations compared to ORN representations (Bhandawat et al., 2007), and 4) a sub-linear transfer function between ORNs and PNs (Bhandawat et al., 2007). Thus, our simulated AL appeared to perform the fundamental computations of real ALs, providing a baseline for assessing the effects of idiosyncratic variation.

Table 2:

We simulated stochastic individuality in the AL circuit in two ways (Figure 4E): 1) glomerular-level variation in PN input-synapse density (reflecting a statistical relationship observed between glomerular volume and synapse density in the hemibrain, Figure 4 – figure supplement 4), and 2) bootstrapping of neuronal compositions within cell types (reflecting variety in developmental program outcomes for ORNs, PNs, etc.). Supplementary Video 4 shows the diverse connectivity matrices attained under these resampling approaches. We simulated odor responses in thousands of ALs made idiosyncratic by these sources of variation, and in each, recorded the firing rates of PNs when stimulated by the 12 odors from our experimental panel (Figure 4F, Figure 4 – figure supplement 1).

To determine which sources of variation produced patterns of PN coding variation consistent with our empirical measurements, we compared principal components of PN responses from real idiosyncratic flies to those of simulated idiosyncratic ALs. Empirical PN responses are strongly correlated at the level of glomeruli (Figure 4G; Figure 1 – figure supplement 8). As a positive control that the model can recapitulate this empirical structure, resampling PN input-synapse density across glomeruli produced PN response correlations strongly organized by glomerulus (Figure 4I). As a negative control, variation in PN responses due solely to poisson timing of ORN input spikes (i.e., absent any circuit idiosyncrasy) was not organized at the glomerular level (Figure 4H). Strikingly, bootstrapping ORN membership yielded a strong glomerular organization in PN responses (Figure 4J). The loadings of the top PCs under ORN bootstrapping are dominated by responses of a single glomerulus to all odors, including DM2 and DC2. This is reminiscent of PC2 of PN calcium responses, with prominent (opposite sign) loadings for DM2 and DC2. Bootstrapping LNs, in contrast, produced much less glomerular organization (Figure 4K), with little resemblance to the loadings of the empirical calcium PCs. The PCA loadings for simulated PN responses under all combinations of cell type bootstrapping and PN input-synapse density resampling are given in Figure 4 – figure supplement 5.

DM2 and DC2 (also DL5) stand out in the PCA loadings under PN input-synapse density resampling and ORN bootstrapping (Figure 4I,J), suggesting that behaviorally-relevant PN coding variation is recapitulated in this modeling framework. To formalize this analysis, for each idiosyncratic AL, we computed a “behavioral preference” by applying the PN PC2 linear model (Figure 1K,M) to simulated PN responses. We then determined how accurately a linear classifier could distinguish OCT- vs MCH-preferring ALs in the space of the first 3 PCs of PN responses (Figure 4 – figure supplement 6). High accuracy was attained under PN input-synapse density resampling and ORN bootstrapping (sources of circuit variation that produced PN response loadings highlighting DM2 and DC2). Thus, developmental variability in ORN populations may drive patterns of PN physiological variation that in turn drive individuality in odor-vs-odor choice behavior.

Discussion

We found an element of the Drosophila olfactory circuit in which patterns of physiological activity predict individual behavioral preferences. This circuit element can be considered a locus of individuality, as it appears to harbor the origins of idiosyncratic preferences among isogenic animals reared in the same environment. Specifically, the difference in the activation of PNs in DC2 and DM2 during odor exposure predicts idiosyncratic OCT-vs-MCH preferences (Figures 1, 2). This circuit element is in the olfactory sensory periphery, suggesting that behavioral idiosyncrasy arises early in the sensorimotor transformation.

Correlating behavior to microscopic circuit features at the individual level is challenging (Koulakov et al., 2005). Measurements of both calcium responses and preference behavior are noisy. Calcium recordings are slow to acquire, making it hard to achieve sample sizes sufficient for machine-learning discovery of correlations with behavior. We conducted three major experiments (paired odor-vs-odor preference and calcium recordings, odor-vs-air preference and calcium recordings, and odor-vs-odor and Brp-Short imaging), each with training and test sets on the scale of 20-60 individuals each. This allowed us to do some limited statistical discovery of correlations, which we restrained by conducting at most five exploratory correlation measurements between circuit and behavioral measures. We were particularly struck by the extent to which PN activity could predict preference between two aversive odors. Importantly, we confirmed this by evaluating the PN calcium-behavior model on a test set of flies measured several weeks after the training flies, finding the same statistically robust trend in both data partitions (training set: R² = 0.15, n = 47, p = 0.0063; testing set: R² = 0.31, n = 22, p = 0.0069; Figure 2 – figure supplement 1).

Previous work has found mammalian peripheral circuit areas are predictive of individual behavior (Britten et al., 1996; Michelson et al., 2017; Newsome et al., 1989; Osborne et al., 2005), but this study is among the first (Linneweber et al., 2020; Mellert et al., 2016; Skutt-Kakaria et al., 2019) to link cellular-level circuit variants and individual behavior in the absence of genetic variation. Another key conclusion is that loci of individuality are likely to vary, even within the sensory periphery, with the specific behavioral paradigm (i.e., odor-vs-odor or odor- vs-air). Our ability to predict behavioral preferences was limited by the repeatability of the behavior itself (Figure 1 – figure supplement 9). Low persistence of odor preference may be attributable to factors like internal states or plasticity. It may be fruitful in future studies to map circuit elements whose activity predicts trial-to-trial behavioral fluctuations within individuals.

Seeking insight into the molecular basis of behaviorally-relevant physiological variation, we imaged Brp in the axon terminals of the ORN-PN synapse, using confocal and expansion microscopy. Brp glomerular density was a significant predictor of individual odor-vs-odor preferences (Figure 3). The strongest predictor of OCT-MCH preference among principal components of Brp-Short density features contrastive loadings between DM2 and other glomeruli, similar to the DM2 - DC2 contrast present in the model that predicts odor preference from PN calcium. This is consistent with the recent finding of a linear relationship between synaptic density and excitatory postsynaptic potentials (Liu et al., 2022) and another study in which idiosyncratic synaptic density in central complex output neurons predicts individual locomotor behavior (Skutt-Kakaria et al., 2019). The predictive relationship between Brp and behavior was weaker than that of PN calcium responses, suggesting there are other determinants, such as other synaptic proteins, neurite morphology, or the influence of idiosyncratic LNs (Chou et al., 2010) modulating the ORN-PN transformation (Nagel et al., 2015).

To integrate our synaptic and physiological results, we implemented a spiking model with 3,062 neurons and synaptic weights drawn directly from the fly connectome (Scheffer et al., 2020) (Figure 4). With light parameter tuning, this model recapitulated canonical AL computations, providing a baseline for assessing the effects of idiosyncratic stochastic variation. The apparent variation in odor responses across simulated individuals (Figure 4F) is less than that seen in the empirical calcium responses (Figure 1H), likely due to 1) biological phenomena missing from the model, 2) the lack of measurement noise, and 3) the fact that our perturbations are applied to the connectome of a single fly. When examining PCA loadings, however, simulating idiosyncratic ALs by varying PN input synapse density or bootstrapping ORNs produced correlated PN responses across odors in DC2 and DM2, matching our experimental results.

These sources of variation specifically implicate the ORN-PN synapse (like our Brp results) as an important substrate for establishing behaviorally-relevant patterns of PN response variation.

The flies used in our experiments were isogenic and reared in standardized laboratory conditions that produce reduced behavioral individuality compared to enriched environments (Akhund-Zade et al., 2019; Körholz et al., 2018; Zocher et al., 2020). Yet, even these conditions yield substantial behavioral individuality. We do not expect variability in the expression of the flies’ transgenes to be a major driver of this individuality, as wildtype flies have a similarly broad distribution of odor preferences (Honegger et al., 2019). The ultimate source of stochasticity in this behavior remains a mystery, with possibilities ranging from thermal fluctuations at the molecular scale to macroscopic, but seemingly irrelevant, variations like the exact fill level of the culture media (Honegger and de Bivort, 2018). Developing nervous systems employ various compensation mechanisms to dampen out the effects of these fluctuations (Marder, 2011; Tobin et al., 2017). Behavioral variation may be beneficial, supporting a bet-hedging strategy (Hopper, 1999) to counter environmental fluctuations (Akhund-Zade et al., 2020; Honegger et al., 2019; Kain et al., 2015; Krams et al., 2021). Empirically, the net effect of dampening and accreted ontological (Gomez-Marin and Ghazanfar, 2019) fluctuations is individuals with diverse behaviors. This process unfolds across all levels of biological regulation. Just as PN response variation appears to be partially rooted in glomerular Brp variation, the latter has its own molecular roots, including, perhaps, stochasticity in gene expression (Li et al., 2017; Raj et al., 2010), itself a predictor of idiosyncratic behavioral biases (Werkhoven et al., 2021). Improved methods to longitudinally assay the fine-scale molecular and anatomical makeup of behaving organisms throughout development and adulthood will be invaluable to further illuminate the mechanistic origins of individuality.

Materials and methods

Data and code availability

All raw data, totaling 600 GB, are available via hard drive from the authors. A smaller (7 GB) repository with partially processed data files and MATLAB/Python scripts sufficient to generate figures and results is available at Zenodo (doi:10.5281/zenodo.8092971).

Fly rearing

Experimental flies were reared in a Drosophila incubator (Percival Scientific DR-36VL) at 22° C, 40% relative humidity, and 12:12h light:dark cycle. Flies were fed cornmeal/dextrose medium, as previously described (Honegger et al., 2019). Mated female flies aged 3 days post-eclosion were used for behavioral persistence experiments. Mated female flies aged 7 to 15 days post-eclosion were used for all paired behavior-calcium imaging and immunohistochemistry experiments.

Fly stocks

The following stocks were obtained from the Bloomington Drosophila Stock Center: P{20XUAS-IVS-GCaMP6m}attP40 (BDSC #42748), w[*]; P{w[+mC]=Or13a-GAL4.F}40.1 (BDSC #9945), w[*]; P{w[+mC]=Or19a-GAL4.F}61.1 (BDSC #9947), w[*]; P{w[+mC]=Or22a-GAL4.7.717}14.2 (BDSC #9951), w[*]; P{w[+mC]=Orco-GAL4.W}11.17; TM2/TM6B, Tb[1] (BDSC #26818). Transgenic lines were outcrossed to the isogenic line isokh11 (Honegger et al., 2019) for at least 5 generations prior to being used in any experiments. GH146-Gal4 was a gift provided by Y. Zhong (Honegger et al., 2019). w; UAS-Brp-Short-mStrawberry; UAS-mCD8-GFP; + was a gift of Timothy Mosca and was not outcrossed to the isokh11 background (Mosca and Luo, 2014).

Odor delivery

Odor delivery during behavioral tracking and neural activity imaging was controlled with isolation valve solenoids (NResearch Inc.) (Honegger et al., 2019). Saturated headspace from 40 ml vials containing 5 ml pure odorant were serially diluted via carbon-filtered air to generate a variably (10-25%) saturated airstream controlled by digital flow controllers (Alicat Scientific) and presented to flies at total flow rates of ∼100 mL/min. Dilution on the order of 10% is typical of other odor tunnel assays, as in Claridge-Chang et al. (2009). To yield the greatest signal of individual odor preference, dilution factors for odorants were adjusted on a week-by-week basis to ensure that the mean preference was approximately 50%. The odor panel used for imaging was comprised of the following odorants: 2-heptanone (CAS #110-43-0, Millipore Sigma), 1-pentanol (CAS #71-41-0, Millipore Sigma), 3-octanol (CAS #589-98-0, Millipore Sigma), hexyl-acetate (CAS #142-92-7, Millipore Sigma), 4-methylcyclohexanol (CAS #589-91-3, Millipore Sigma), pentyl acetate (CAS #628-63-7, Millipore Sigma), 1-butanol (CAS #71-36-3, Millipore Sigma), ethyl lactate (CAS #97-64-3, Millipore Sigma), geranyl acetate (CAS #105-87-3, Millipore Sigma), 1-hexanol (CAS #111-27-34, Millipore Sigma), citronella java essential oil ( 191112, Aura Cacia), and 200 proof ethanol (V1001, Decon Labs).

Odor preference behavior

Odor preference was measured at 25°C and 20% relative humidity. As previously described (Honegger et al., 2019), individual flies confined to custom-fabricated tunnels were illuminated with infrared light and behavior was recorded with a digital camera (Basler) and zoom lens (Pentax). The odor choice tunnels were 50 mm long, 5 mm wide, and 1.3 mm tall. Custom real-time tracking software written in MATLAB was used to track centroid, velocity, and principal body axis angle throughout the behavioral experiment, as previously described (Honegger et al., 2019). After a 3-minute acclimation period, odorants were delivered to either end of the tunnel array for 3 minutes. Odor preference score was calculated as the fraction of time spent in the reference side of the tunnel during odor-on period minus the time spent in the reference side of the tunnel during the pre-odor acclimation period.

Behavioral preference persistence measurements

After measuring odor preference, flies were stored in individual housing fly plates (modified 96-well plates; FlySorter, LLC) on standard food, temperature, humidity, and lighting conditions.

Odor preference of the same individuals was measured 3 and/or 24 hours later. In some cases, fly tunnel position was randomized between measurements. We observed that randomization had little effect on preference persistence.

Calcium imaging

Flies expressing GCaMP6m in defined neural subpopulations were imaged using a custom-built two-photon microscope and ultrafast Ti:Sapphire laser (Spectra-Physics Mai Tai) tuned to 930 nm, at a power of 20 mW out of the objective (Olympus XLUMPlanFL N 20x/1.00 W). For paired behavior and imaging experiments, the time elapsed between behavior measurement and imaging ranged from 15 minutes to 3 hours. Flies were anesthetized on ice and immobilized in an aluminum sheet with a female-fly-sized hole cut in it. The head cuticle between the antennae and ocelli was removed along with the tracheae to expose the ALs from the dorsal side. Volume scanning was performed using a piezoelectric objective mount (Physik Instrumente). ScanImage 2013 software (Vidrio Technologies) was used to coordinate galvanometer laser scanning and image acquisition. Custom Matlab (Mathworks) scripts were used to coordinate image acquisition and control odor delivery. 256 by 192 (x-y) pixel 16-bit tiff images were recorded.

The piezo travel distance was adjusted between 70 and 90 μm so as to cover most of the AL. The number of z-sections in a given odor panel delivery varied between 7 and 12 yielding a volume acquisition rate of 0.833 Hz. Odor delivery occurred from 6-9.6s of each recording.

Each fly experienced up to four deliveries of the odor panel. The antennal lobe being recorded (left or right) was alternated after each successful completion of an odor panel. Odors were delivered in randomized order. In cases where baseline fluorescence was very weak or no obvious odor responses were visible, not all four panels were delivered.

Glomerulus segmentation and labeling

Glomerular segmentation masks were extracted from raw image stacks using a k-means clustering algorithm based on time-varying voxel fluorescence intensities, as previously described (Honegger et al., 2019). Each image stack, corresponding to a single odor panel delivery, was processed individually. Time-varying voxel fluorescence values for each odor delivery were concatenated to yield a voxel-by-time matrix consisting of each voxel’s recorded value during the course of all 13 odor deliveries of the odor panel. After z-scoring, principal component analysis was performed on this matrix and 75% of the variance was retained. Next, k-means (k = 80, 50 replicates with random starting seeds) was performed to produce 50 distinct voxel cluster assignment maps which we next used to calculate a consensus map. This approach was more accurate than clustering based on a single k-means seed.

Of the 50 generated voxel cluster assignment maps, the top 5 were selected by choosing those maps with the lowest average within-cluster sum of distances, selecting for compact glomeruli. The remaining maps were discarded. Next, all isolated voxel islands in each of the top 5 maps were identified and pruned based on size (minimum size = 100 voxels, maximum size = 10000 voxels). Finally, consensus clusters were calculated by finding voxel islands with significant overlap across all 5 of the pruned maps. Voxels which fell within a given cluster across all 5 pruned maps were added to the consensus cluster. This process was repeated for all clusters until the single consensus cluster map was complete. In some cases we found by manual inspection that some individual glomeruli were clearly split into two discrete clusters. These splits were remedied by automatically merging all consensus clusters whose centroids were separated by a physical distance of less than 30 voxels and whose peak odor response Spearman correlation was greater than 0.8. Finally, glomeruli were manually labeled based on anatomical position, morphology, and size (Grabe et al., 2015). We focused our analysis on 5 glomeruli (DM1, DM2, DM3, DL5, and DC2), which were the only glomeruli that could be observed in all paired behavior-calcium datasets. However, not all 5 glomeruli were identified in all recordings (Figure 1 – figure supplement 3). Missing glomerular data was later mean-imputed. Using alternating least squares to impute data (running the pca function with option ‘als’ to infill missing values 1,000 times and taking the mean infilled matrix – see Figure 1 – figure supplement 5 of Werkhoven et al., 2021) had negligible effect on the fitted slope and predictive capacity of the PN PC2 OCT-MCH model compared to mean-infilling.

Calcium image data analysis

All data was processed and analyzed in MATLAB 2018a (Mathworks). Calcium responses for each voxel were calculated as Δf/f = [f(t) - F]/F, where f(t) and F are the instantaneous and average fluorescence, respectively. Each glomerulus’ time-dependent calcium response was calculated as the mean Δf/f across all voxels falling within the glomerulus’ automatically-generated segmentation mask during a single volume acquisition. Time-varying odor responses were normalized to baseline by subtracting the median of pre-odor Δf/f from each trace. Peak odor response was calculated as the maximum fluorescence signal from 7.2s to 10.8s (images 6 through 9) of the recording.

To compute principal components of calcium dynamics, each fly’s complement of odor panel responses (a 5 glomeruli by 13 odors = 65-dimensional vector) was concatenated. Missing glomerulus-odor response values were filled in with the mean glomerulus-odor pair across all fly recordings for which the data was not missing. After infilling, principal component analysis was carried out with individual odor panel deliveries as observations and glomerulus-odor responses pairs as features.

Inter- and intra-fly distances (Figure 1J) were calculated using the projections of each fly’s glomerulus-odor responses onto all principal components. For each fly, the average Euclidean distance between response projections 1) among left lobe trials, 2) among right lobe trials, and 3) between left and right lobe trials were averaged together to get a single within-fly distance. Intra-fly distances were computed in a similar fashion (for each fly, taking the average distance of its response projections to those of other flies using only left lobe trials / only right lobe trials / between left-right trials, then averaging these three values to get a single across-fly distance).

In a subset of experiments in which we imaged calcium activity, some solenoids failed to open, resulting in the failure of odor delivery in a small number of trials. In these cases, we identified trials with valve failures by manually recognizing that glomeruli failed to respond during the nominal odor period. These trials were treated as missing data and infilled, as described above. Fewer than ∼10% of flies and 5% of odor trials were affected.

For all predictive models constructed, the average principal component score or glomerulus-odor Δf/f response across trials was used per individual; that is, each fly contributed one data point to the relevant model. Linear models were constructed from behavior scores and the relevant predictor (principal component, average Δf/f across dimensions, specific glomerulus measurements) as described in the text and Tables 1-2. All reported linear model p-values are nominal, i.e., unadjusted for multiple hypothesis comparisons. 95% confidence intervals around model regression lines were estimated as +/- 2 standard deviations of the value of the regression line at each x-position across 2000 bootstrap replicates (resampling flies). To predict behavior as a function of time during odor delivery, we analyzed data as described above, but considered only Δf/f at each single time point (Figure 2 – figure supplement 2A-C), rather than averaging during the peak response interval.

To decode individual identity from neural responses, we first performed PCA on individual odor panel peak responses. We retained principal component scores constituting specified fractions of variance (Figure 1 – figure supplement 6A) and trained a linear logistic classifier to predict individual identity from single odor panel deliveries.

To decode odor identity from neural responses, each of the 5 recorded glomeruli were used as features, and the calcium response of each glomerulus to a specific odor at a specified time point were used as observations (PNs, n = 5317 odor deliveries; ORNs, n = 2704 odor deliveries). A linear logistic classifier was trained to predict the known odor identity using 2-fold cross-validation. That is, a model was trained on half the data and evaluated on the remaining half, and then this process was repeated with the train and test half reversed. The decoding accuracy was quantified as the fraction of odor deliveries in which the predicted odor was correct.

Inference of correlation between latent calcium and behavior states

We performed a simulation-based analysis to infer the strength of the correlation between latent calcium (Brp) and behavior states, given the R² of a given linear model. Figure 1 – figure supplement 9 is a schematic of a possible data generation process that underlies our observed data. We assume that the “true” behavioral and calcium values of the animal are captured by unobserved latent states X_c and X_b, respectively, such that the correlation between X_c and X_b is the biological signal captured by the model, having adjusted for the noise associated with actually measuring behavior and calcium (⍴_signal ). Our calcium and odor preference scores are subject to measurement error and temporal instability (behavior and neural activity were measured 1-3 hours apart). These effects are both noise with respect to estimating the linear relationship between calcium and behavior. Their magnitude can be estimated using the empirical repeatability of behavior and calcium experiments respectively. Thus, our overall approach was to assume true latent behavior and calcium signals that are correlated at the level ⍴_signal, add noise to them commensurate with the repeatability of these measures to simulate measured behavior and calcium, and record the simulated empirical R² between these measured signals. This was done many times to estimate distributions of empirical R² given ⍴_signal. These distributions could finally be used in the inverse direction to infer ⍴_signal given the actual model R² values computed in our study.

Specifically, we simulated X_c as a set of N standard normal variables (N equalling the number of flies used to compute a correlation between predicted and measured preference) and generated X_b = ⍴_signal X_c + (1-⍴_signal²Z)^½, where Z is a set of N standard normal variables uncorrelated with X_c, a procedure that ensures that corr(X_c, X_b) = ⍴_signal. Next, we simulated observed calcium readouts X_c’ and X_c”, such that corr(X_c, X_c’) = corr(X_c, X_c”) = r_c. Similarly, we simulated noisy observed behavioral assay readouts X_b’ and X_b”, such that corr(X_b, X_b’) = corr(X_b, X_b”) = r_b. The values of r_c and r_b were drawn from the empirical repeatability of calcium (R_c,c²) and behavior (R_b,b²) respectively as follows. Since calcium is a multidimensional measure, and our calcium model predictors are based on principal components of glomerulus-odor responses, we used variance explained along the PCs to calculate a single value for the calcium repeatability R_c,c². We compared the eigenvalues of the real calcium PCA to those of shuffled calcium data (shuffling glomerulus/odor responses for each individual fly), computing R_c,c² by summing the variance explained along the PCs of the calcium data up until the component-wise variance for the calcium data fell below that of the shuffled data, a similar approach as done in Berman et al., 2014 and Werkhoven et al., 2021. That is, we determined which empirical PCs had more variance than their corresponding rank-matched PC in shuffled data, interpreted the remaining PCs as harboring the noise of the experiment, and totaled the variance explained of the non-noise PCs as our measure of the repeatability of the measurement as a whole. R_c,c² was calculated to be 0.77 for the full PN calcium data.

To incorporate uncertainty in calcium-calcium repeatability, we utilized bootstrapping. We resampled the calcium data associated with individual flies 10,000 times, performed PCA and computed R_c,c² for each resampled dataset, then set r_c = (R_c,c²)^1/4 to ensure corr(X_c’, X_c”)² = R_c,c². For behavior-behavior uncertainty, we set r_b from the repeatability across odor preference trials in the same flies measured 3h apart (R_b,b² = 0.12 for OCT vs MCH, Figure 1 – figure supplement 1D using the full dataset of flies). We also resampled the flies 10,000 times, computed R_b,b² for each resampled dataset, and set r_b = (R_b,b²)^1/4 to ensure corr(X_b’, X_b”)² = R_b,b².

We varied ⍴_signal from 0 to 1 in increments of 0.01, and for each ⍴_signal and bootstrap iteration we simulated a set of N X_c, and generated X_b , X_c’, X_c”, X_b’, and X_b”, then we computed a simulated observed calcium-behavior relationship strength R_c,b² = corr(X_c’, X_b’)². We repeated this simulation 10,000 times for each ⍴_signal and plotted the resultant relationship between ⍴_signal against R_c,b² (percentiles of R_c,b² are displayed in Figure 1 – figure supplement 9B). Then, we inferred ⍴_signal by first drawing bootstrapped samples of calcium-behavior R², then adding together the marginal distributions of ⍴_signal for each calcium-behavior R². We report the median ⍴_signal and 90% confidence interval as estimated by the 5th-95th quantiles.

The procedure outlined above was done analogously for models using Brp-short relative fluorescence intensity, performing the PCA-based calcium response repeatability step with PCA on the multidimensional Brp-short relative fluorescence intensity (which yielded R_brp,brp² = 0.78).

DoOR data

DoOR data for the glomeruli and odors relevant to our study was downloaded from http://neuro.uni-konstanz.de/DoOR/default.html (Münch and Galizia, 2016).

Yoked odor experience experiments

We selected six flies for which both odor preference and neural activity were recorded to serve as the basis for imposed odor experiences for yoked control flies. The experimental flies were chosen to represent a diversity of preference scores. Each experimental fly’s odor experience was binned into discrete odor bouts to represent experience of either MCH or OCT based on its location in the tunnel as a function of time (Figure 2D). Odor bouts lasting less than 100 ms were omitted due to limitations on odor-switching capabilities of the odor delivery apparatus. To deliver a given experimental fly’s odor experience to yoked controls, we set both odor streams (on either end of the tunnel apparatus) to deliver the same odor experienced by the experimental fly at that moment during the odor-on period. No odor was delivered to yoked controls during time points in which the experimental fly resided in the tunnel choice zone (central 5 mm). See Figure 2D for an example pair of experimental fly and yoked control behavior and odor experience.

Immunohistochemistry

After measuring odor preference behavior, 7-15 day-old flies were anesthetized on ice and brains were dissected in phosphate buffered saline (PBS). Dissection and immunohistochemistry were carried out as previously reported (Wu and Luo, 2006). The experimenter was blind to the behavioral scores of all individuals throughout dissection, imaging, and analysis. Individual identities were maintained by fixing, washing, and staining each brain in an individual 0.2 mL PCR tube using fluid volumes of 100 uL per brain (Fisher Scientific). Primary incubation solution contained mouse anti-nc82 (1:40, DSHB), chicken anti-GFP (1:1000, Aves Labs), rabbit anti-mStrawberry (1:1000, biorbyt), and 5% normal goat serum (NGS, Invitrogen) in PBT (0.5% Triton X-100 in PBS). Secondary incubation solution contained Atto 647N-conjugated goat anti-mouse (1:250, Millipore Sigma), Alexa Fluor 568-conjugated goat anti-rabbit (1:250), Alexa Fluor 488-conjugated goat anti-chicken (1:250, ThermoFisher), and 5% NGS in PBT. Primary and secondary incubation times were 2 and 3 overnights, respectively, at 4° C. Stained samples were mounted and cleared in Vectashield (H-1000, Vector Laboratories) between two coverslips (12-568B, Fisher Scientific). Two reinforcement labels (5720, Avery) were stacked to create a 0.15 mm spacer.

Expansion microscopy

Immunohistochemistry for expansion microscopy was carried out as described above, with the exception that antibody concentrations were modified as follows: mouse anti-nc82 (1:40), chicken anti-GFP (1:200), rabbit anti-mStrawberry (1:200), Atto 647N-conjugated goat anti-mouse (1:100), Alexa Fluor 568-conjugated goat anti-rabbit (1:100), Alexa Fluor 488-conjugated goat anti-chicken (1:100). Expansion of stained samples was performed as previously described (Asano et al., 2018; Gao et al., 2019). Expanded samples were mounted in coverslip-bottom petri dishes (MatTek Corporation) and anchored by treating the coverslip with poly-l-lysine solution (Millipore Sigma) as previously described (Asano et al., 2018).

Confocal imaging

All confocal imaging was carried out at the Harvard Center for Biological Imaging. Unexpanded samples were imaged on an LSM700 (Zeiss) inverted confocal microscope equipped with a 40x oil-immersion objective (1.3 NA, EC Plan Neofluar, Zeiss). Expanded samples were imaged on an LSM880 (Zeiss) inverted confocal microscope equipped with a 40x water-immersion objective (1.1 NA, LD C-Apochromat, Zeiss). Acquisition of Z-stacks was automated with Zen Black software (Zeiss).

Standard confocal image analysis

We used custom semi-automated code to generate glomerular segmentation masks from confocal z-stacks of unexpanded Orco>Brp-Short brains. Using Matlab, each image channel was median filtered (σ_x, σ_y, σ_z = 11, 11, 1 pixels) and downsampled in x and y by a factor of 11. Next, an ORN mask was generated by multiplying and thresholding the Orco>mCD8 and Orco>Brp-Short channels. Next, a locally normalized nc82 and Orco>mCD8 image stack were multiplied and thresholded, and the ORN mask was applied to remove background and other undesired brain structures. This pipeline resulted in a binary image stack which maximized the contrast of the glomerular structure of the antennal lobe. We then applied a binary distance transform and watershed transform to generate discrete subregions which aimed to represent segmentation masks for each glomerulus tagged by Orco-Gal4.

However, this procedure generally resulted in some degree of under-segmentation; that is, some glomerular segmentation masks were merged. To split each merged segmentation mask, we convolved a ball (whose radius was proportional to the cube root of the volume of the segmentation mask in question) across the mask and thresholded the resulting image. The rationale of this procedure was that 2 merged glomeruli would exhibit a mask shape resembling two touching spheres, and convolving a similarly-sized sphere across this volume followed by thresholding would split the merged object. After ball convolution, we repeated the distance and watershed transform to once more generate discrete subregions representing glomerular segmentation masks. This second watershed step generally resulted in over-segmentation; that is, by visual inspection it was apparent that many glomeruli were split into multiple subregions.

Therefore, we finally manually agglomerated the over-segmented subregions to generate single segmentation masks for each glomerulus of interest. We used a published atlas to aid manual identification of glomeruli (Grabe et al., 2015). The total Brp-Short fluorescence signal within each glomerulus was determined and divided by the volume of the glomerulus’ segmentation mask to calculate Brp-Short density values.

Expansion microscopy image analysis

The spots function in Imaris 9.0 (Bitplane) was used to identify individual Brp-Short puncta in expanded sample image stacks of Or13a>Brp-Short samples (Mosca and Luo, 2014). The spot size was set to 0.5 um, background subtraction and region-growing were enabled, and the default spot quality threshold was used for each image stack. Identified spots were used to mask the Brp-Short channel and the resultant image was saved as a new stack. In MATLAB, a glomerular mask was generated by smoothing (σ_x, σ_y, σ_z = 40, 40, 8 pixels) and thresholding (92.5th percentile) the raw Brp-Short image stack. The mask was then applied to the spot image stack to remove background spots. Finally, the masked spot image stack was binarized and spot number and properties were quantified.

Antennal Lobe modeling

We constructed a model of the antennal lobe to test the effect of circuit variation on PN activity variation across individuals. Our general approach to producing realistic circuit activity with the AL model was 1) using experimentally-measured parameters whenever possible (principally the connectome wiring diagram and biophysical parameters measured electrophysiologically), 2) associating free parameters only with biologically plausible categories of elements, while minimizing their number, and 3) tuning the model using those free parameters so that it reproduced high-level patterns of activity considered in the field to represent the canonical operations of the AL. Simulations were run in Python (version 3.6) (van Rossum and Drake, 2011), and model outputs were analyzed using Jupyter notebooks (Kluyver et al., 2016) and Python and Matlab scripts.

AL model neurons

Release 1.2 of the hemibrain connectomics dataset (Scheffer et al., 2020) was used to set the connections in the model. Hemibrain body IDs for ORNs, LNs, and PNs were obtained via the lists of neurons supplied in the supplementary tables in Schlegel et al., 2020. ORNs and PNs of non-olfactory glomeruli (VP1d, VP1l, VP1m, VP2, VP3, VP4, VP5) were ignored, leaving 51 glomeruli. Synaptic connections between the remaining 2574 ORNs, 197 LNs, 166 mPNs, and 130 uPNs were queried from the hemibrain API. All ORNs were assigned to be excitatory (Wilson, 2013). Polarities were assigned to PNs based on the neurotransmitter assignments in Bates et al., 2020. mPNs without neurotransmitter information were randomly assigned an excitatory polarity with probability equal to the fraction of neurotransmitter-identified mPNs that are cholinergic; the same process was performed for uPNs. After confirming that the model’s output was qualitatively robust to which mPNs and uPNs were randomly chosen, this random assignment was performed once and then frozen for subsequent analyses.

Of the 197 LNs, we assigned 31 to be excitatory, based on the estimated 1:5.4 ratio of eLNs to iLNs in the AL (Tsai et al., 2018). To account for observations that eLNs broadly innervate the AL (Shang et al., 2007), all LNs were ranked by the number of innervated glomeruli, and the 31 eLNs were chosen uniformly at random from the top 50% of LNs in the list. This produced a distribution of glomerular innervations in eLNs qualitatively similar to that of krasavietz LNs in Supplementary Figure 6 of Chou et al., 2010.

Voltage model

We used a single-compartment leaky-integrate-and-fire voltage model for all neurons as in Kakaria and de Bivort, 2017, in which each neuron had a voltage V_i(t) and current I_i(t). When the voltage of neuron i was beneath its threshold V_{i, thr}, the following dynamics were obeyed:

Each neuron i had electrical properties: membrane capacitance C_i, resistance R_i, and resting membrane potential V_i,0 with values from electrophysiology measurements (Table 2).

When the voltage of a neuron exceeded the threshold V_{i, thr}, a templated action potential was filled into its voltage time trace, and a templated postsynaptic current was added to all downstream neurons, following the definitions in Kakaria and de Bivort, 2017.

Odor stimuli were simulated by triggering ORNs to spike at frequencies matching known olfactory receptor responses to the desired odor. The timing of odor-evoked spikes was given by a Poisson process, with firing rate FR for ORNs of a given glomerulus governed by:

FR_max, the maximum ORN firing rate, was set to 400 Hz. D_{glom, odor} is a value between 0 and 1 from the DoOR database, representing the response of an odorant receptor/glomerulus to an odor, estimated from electrophysiology and/or fluorescence data (Münch and Galizia, 2016). ORNs display adaptation to odor stimuli (Wilson, 2013), captured by the final term with timescale t_a = 110 ms to 75% of the initial value, as done in Kao and Lo, 2020. Thus, the functional maximum firing rate of an ORN was 75% of 400 Hz = 300 Hz, matching the highest ORN firing rates observed experimentally (Hallem et al., 2004). After determining the times of ORN spikes according to this firing-rate rule, spikes were induced by the addition of 10⁶ picoamps in a single time step. This reliably triggered an action potential in the ORN, regardless of currents from other neurons. In the absence of odors, spike times for ORNs were drawn by a Poisson process at 10 Hz, to match reported spontaneous firing rates (de Bruyne et al., 2001).

For odor-glomeruli combinations with missing DoOR values (40% of the dataset), we performed imputation via alternating least squares (using the pca function with option ‘als’ to infill missing values (MATLAB documentation) on the odor x glomerulus matrix 1000 times and taking the mean infilled matrix, which provides a closer match to ground truth missing values than a single run of ALS (Figure 1 – figure supplement 5 of Werkhoven et al., 2021).

A neuron j presynaptic to i supplies its current I_j(t) scaled by the synapse strength W_ji, the number of synapses in the hemibrain dataset from neuron j to i. Rows in W corresponding to neurons with inhibitory polarity (i.e. GABAergic PNs or LNs) were set negative. Finally, post-synaptic neurons (columns of the connectivity matrix) have a class-specific multiplier a_i, a hand-tuned value, described below.

AL model tuning

Class-specific multiplier current multipliers (a_i) were tuned using the panel of 18 odors from Bhandawat et al., 2007 (our source for several experimental observations of high-level AL function): benzaldehyde, butyric acid, 2,3-butanedione, 1-butanol, cyclohexanone, Z3-hexenol, ethyl butyrate, ethyl acetate, geranyl acetate, isopentyl acetate, isoamyl acetate, 4-methylphenol, methyl salicylate, 3-methylthio-1-propanol, octanal, 2-octanone, pentyl acetate, E2-hexenal, trans-2-hexenal, gamma-valerolactone. Odors were “administered” for 400 ms each, with 300 ms odor-free pauses between odor stimuli.

The high-level functions of the AL that represent a baseline, working condition were: (1) firing rates for ORNs, LNs, and PNs matching the literature (listed in Table 2 and see (Bhandawat et al., 2007; Dubin and Harris, 1997; Jeanne and Wilson, 2015; Seki et al., 2010), (2) a more uniform distribution of PN firing rates during odor stimuli compared to ORN firing rates, (3) greater separation of representations of odors in PN-coding space than in ORN-coding space, and (4) a sublinear transfer function between ORN firing rates and PN firing rates. Features (2) - (4) relate to the role of the AL in enhancing the separability of similar odors (Bhandawat et al., 2007).

To find a parameterization with those functions, we tuned the values of a_i as scalar multipliers on ORN, eLN, iLN, and PN columns of the hemibrain connectivity matrix. Thus, these values represent cell type-specific sensitivities to presynaptic currents, which may be justified by the fact that ORNs/LNs/PNs are genetically distinct cell populations (McLaughlin et al., 2021; Xie et al., 2021). A grid search of the four class-wise sensitivity parameters produced a configuration that reasonably satisfied the above criteria (Figure 4 – figure supplement 2). In this configuration, the ORN columns of the hemibrain connectivity matrix are scaled by 0.1, eLNs by 0.04, iLNs by 0.02, and PNs by 0.4. The relatively large multiplier on PNs is potentially consistent with the fact that PNs are sensitive to small differences between weak ORN inputs (Bhandawat et al., 2007). Model outputs were robust over several different sets of a_i, provided iLN sensitivity ≃ eLN < ORN < PN.

We analyzed the sensitivity of the model’s parameters around their baseline values of a_ORN, a_eLN, a_iLN, a_PN = (0.1, 0.04, 0.02, 0.4). Each parameter was independently scaled up to 4x or 1/4x of its baseline value (Figure 4 – figure supplement 3), and the PN firing rates recorded. Separately, multiple-parameter manipulations were performed by multiplying each parameter by a random log-Normal value with mean 1 and +/-1 standard deviation corresponding to a 2x or 0.5x scaling on each parameter. Mean PN-odor responses were calculated for all manipulated runs and compared to the mean PN-odor responses for the baseline configuration. A manipulation effect size was calculated by cohen’s d ((mean manipulated response - mean baseline response)/(pooled standard deviation)). None of these manipulations reached effect size magnitudes larger than 0.9 (which can be roughly interpreted as the number of standard deviations in the baseline PN responses away from the mean baseline PN response), which signaled that the model was robust to the sensitivity parameters in this range. The most sensitive parameter was, unsurprisingly, a_PN.

Notable ways in which the model behavior deviates from experimental recordings (and thus caveats on the interpretation of the model) include: 1) Model LNs appear to have more heterogeneous firing rates than real LNs, with many LNs inactive for this panel of odor stimuli. This likely reflects a lack of plastic/homeostatic mechanisms in the model to regularize LN firing rates given their variable synaptic connectivity (Chou et al., 2010). 2) Some PNs had off-odor rates that are high compared to real PNs, resulting in a distribution of ON-OFF responses that had a lower limit than in real recordings. Qualitatively close matches were achieved between the model and experimental data in the distributions of odor representations in ORN vs PN spaces and the non-linearity of the ORN-PN transfer function.

AL model circuit variation generation

We generated AL circuit variability in two ways: cell-type bootstrapping, and synapse density resampling. These methods assume that the distribution of circuit configurations across individual ALs can be generated by resampling circuit components within a single individual’s AL (neurons and glomerular synaptic densities, respectively, from the hemibrain EM volume).

To test the effect of variation in the developmental complement of neurons of particular types, we bootstrapped populations of interest from the list of hemibrain neurons. Resampling with replacement of ORNs was performed glomerulus-by-glomerulus, i.e., separately among each pool of ORNs expressing a particular Odorant receptor gene. The same was done for PNs. For LNs, all 197 LNs were treated as a single pool; there was no finer operation based on LN subtypes or glomerular innervations. This choice reflects the high developmental variability of LNs (Chou et al., 2010). The number of synapses between a pair of bootstrapped neurons was equal to the synapse count between those neurons in the hemibrain connectivity matrix.

In some glomeruli, bootstrapping PNs produced unreasonably high variance in the total PN synapse count. For instance, DP1m, DC4, and DM3 each harbor PNs that differ in total synapse count by a factor of ∼10. Since these glomeruli have between two to three PNs each, in a sizable proportion of bootstrap samples, all-highly connected (or all-lowly) connected PNs are chosen in such glomeruli. To remedy this biologically unrealistic outcome, we examined the relationship between total input PN synapses within a glomerulus and glomerular volume (Figure 4 – figure supplement 4). In the “synapse density resampling” method, we required that the number of PN input synapses within a glomerulus reflect a draw from the empirical relationship between total input PN synapses and glomerular volume as present in the hemibrain data set. This was achieved by, for each glomerulus, sampling from the following distribution that depends on glomerular volume, then multiplying the number of PN input synapses by a scalar to match that sampled value:

Here S_g is the PN input synapse count for glomerulus g, V_g is the volume of glomerulus g (in cubic microns), ε is a Gaussian noise variable with standard deviation σ, and a, d are the scaling factor and exponent of the volume term, respectively. The values of these parameters (a = 8.98, d = 0.73, σ = 0.38) were fit using maximum likelihood.

Quantification and statistical analysis

All fly behavior and calcium data was processed and analyzed in MATLAB 2018a (Mathworks). AL simulations were run in Python (version 3.6) (van Rossum and Drake, 2011), and model outputs were analyzed using Jupyter notebooks (Kluyver et al., 2016) and Python scripts. We performed a power analysis prior to the study to determine that recording calcium activity in 20-40 flies would be sufficient to identify moderate calcium-behavior correlations. Sample sizes for expansion microscopy were smaller, as the experimental procedure was more involved – therefore, we did not conduct a formal statistical analysis. Linear models were fit using the fitlm MATLAB function (https://www.mathworks.com/help/stats/fitlm.html); coefficients and p-values of models between measured preferences and predicted preferences are listed in Table 1. 95% confidence intervals around model regression lines were estimated as +/- 2 standard deviations of the value of the regression line at each x-position across 2000 bootstrap replicates (resampling flies). Boxplots depict the median value (points), interquartile range (boxes), and range of the data (whiskers).

Acknowledgements

We thank Asa Barth-Maron and Rachel Wilson for discussions helpful to the AL modeling, and Katrin Vogt for help revising the manuscript. Ed Soucy and Brett Graham of the Center for Brain Science Neuroengineering Core helped maintain the olfactometer and microscope. D.L. was supported by the NSF-Simons Center for Mathematical and Statistical Analysis of Biology at Harvard, award number #1764269 and the Harvard Quantitative Biology Initiative. B.L.d.B. was supported by a Klingenstein-Simons Fellowship Award, a Smith Family Odyssey Award, a Harvard/MIT Basic Neuroscience Grant, National Science Foundation grant no. IOS-1557913, and NIH/NINDS grant no. 1R01NS121874-01. E.B. was supported by a Harvard/MIT Basic Neuroscience Grant, Lisa Yang, John Doerr, and NIH grant no. 1R01EB024261.

Additional information

Author contributions

M.C. conducted behavior experiments with assistance from M.S., conducted confocal microscopy, and conducted expansion microscopy with contributions from R. G. and E.B. D.L. implemented the computational AL model. B.L.d.B. supervised the project.

Declaration of interests

E.B. is a co-founder of a company that aims to commercialize expansion microscopy for medical purposes. R.G. and E.B. are co-inventors on multiple patents related to expansion microscopy. The authors declare no other competing interests.

Figure supplements

Figure 1 – figure supplement 1.

Figure 1 – figure supplement 2.

Figure 1 – figure supplement 3.

Figure 1 – figure supplement 4.

Figure 1 – figure supplement 5.

Figure 1 – figure supplement 6.

Figure 1 – figure supplement 7.

Figure 1 – figure supplement 8.

Figure 1 – figure supplement 9.

Figure 1 – figure supplement 10.

Figure 2 – figure supplement 1.

Figure 2 – figure supplement 2.

Figure 3 – figure supplement 1.

Figure 3 – figure supplement 2.

Figure 4 – figure supplement 1.

Figure 4 – figure supplement 2.

Figure 4 – figure supplement 3.

Figure 4 – figure supplement 4.

Figure 4 – figure supplement 5.

Figure 4 – figure supplement 6.