Introduction

Understanding the intricate relationship between brain function and energy metabolism is a fundamental challenge in neuroscience. Cellular metabolism relies on extremely conserved biochemical processes that sequentially catabolize primary energy substrates into usable energy in the form of ATP. These ubiquitous processes are especially active in the brain, which consumes large amounts of energy relative to its weight (Rolfe and Brown, 1997), mostly in the form of glucose (Sokoloff, 1999). Neurons indeed have high and dynamic energy needs, due in particular to the maintenance of ion gradients and the propagation of action potentials (Attwell and Laughlin, 2001).

In vertebrates, astrocytes are major players in the achievement of this spatiotemporally fluctuating energy allocation. Astrocytes orchestrate the neurovascular coupling (Bélanger et al., 2011) but also store energy-carrying molecules and fuel them to neurons, thanks to their complementary metabolic profile. While neurons predominantly generate ATP through oxidative phosphorylation, astrocytes display a high level of glycolysis (Itoh et al., 2003; Boumezbeur et al., 2010), resulting, for a significant fraction, in lactate production (Serres et al., 2005; Bouzier-Sore et al., 2006). Numerous studies tend to demonstrate the net transfer of energy from astrocytes to neurons through a so-called lactate shuttle (Pellerin and Magistretti, 1994, 2012; Bélanger et al., 2011). Despite some differences, accumulating evidence shows that insect glia shares many properties with its vertebrate counterparts (Nagai et al., 2021; De Backer and Grunwald Kadow, 2022), including their highly glycolytic profile and their ability to transfer glycolysis-derived metabolites to neurons (Volkenhoff et al., 2015; Rittschof and Schirmeier, 2018).

Additionally to the necessity of adapting energy supply to meet neuronal consumption, studies in diverse species, including insects (Chandrasekaran et al., 2015; Plaçais et al., 2017) and rodents (Shimizu et al., 2007; Parsons and Hirasawa, 2010; Suzuki et al., 2011), have shown that energy fluxes modulate higher brain functions. These findings support a shift away from a purely neurocentric perspective, suggesting that energy is not only delivered to neurons on demand but is instead regulated by an integrated neuron-glia network, actively shaping cerebral activity. Memory formation, in particular, has already been shown to be closely associated with neuronal energy regulation. Several studies have for example reported the importance of astrocytic metabolism during memory formation in vertebrates (Gibbs et al., 2006; Newman et al., 2011; Suzuki et al., 2011; Gao et al., 2016).

With its well-characterized genetics, neuroanatomy and memory capacities, Drosophila serves as an ideal model organism for investigating the neural basis of memory. Aversive olfactory memory can be induced in Drosophila through classical Pavlovian conditioning, by pairing an odor with electric shocks (Tully and Quinn, 1985; Modi et al., 2020). The neurons encoding these memories are located in the mushroom bodies (MBs), which are considered as the main brain center for associative memory in insects (Heisenberg, 2003; Modi et al., 2020). The Drosophila MBs are paired structures including around 2,000 intrinsic neurons, named Kenyon cells (KCs), in each hemisphere (Aso et al., 2009, 2014; Li et al., 2020). These neurons receive dendritic inputs from the antennal lobes in the MB calyces (Keller and Vosshall, 2003), located in the posterior part of the brain. Their axons form a bundle known as the peduncle, which extends to the anterior part of the brain before branching into horizontal and vertical lobes. KCs are categorized into three main subtypes based on their axonal morphology and genetic expression profile: α/β, α’/β’ and γ neurons (Crittenden et al., 1998; Lee et al., 1999). All three populations are considered to be involved at multiple stages of memory processing (Krashes et al., 2007), but converging evidence suggests that γ KCs are required for encoding and retrieving immediate memory (Qin et al., 2012; Bouzaiane et al., 2015; Owald et al., 2015; Perisse et al., 2016), while α/β neurons underlie a middle-term memory that develops more slowly over several hours (Blum et al., 2009; Séjourné et al., 2011; Bouzaiane et al., 2015).

Recent studies have underscored the importance of energy regulation in the MBs during memory formation. Increased mitochondrial energy flux in MB neurons (Plaçais et al., 2017) with precise spatiotemporal dynamics (Comyn et al., 2024) and fueled by the transfer of glycolysis-derived alanine from cortex glia or astrocytes to neurons (Rabah et al., 2023) indeed play critical roles in the formation of middle-term and long-term memory. These findings position the Drosophila MBs as a compelling model for exploring how neuronal metabolism and cognitive processes intersect. In these works, metabolic dynamics in the MBs have been tracked using genetically encoded fluorescent sensors tracking specific metabolites.

The growing recognition of the role played by energy metabolism in brain function highlights the need for universal methods to monitor cellular metabolic profiles in vivo across species, with high spatial and temporal resolution. Cellular ATP production through catabolic pathways relies on redox reactions involving carbon-based energy substrates and redox couples, notably the coenzyme nicotinamide adenine dinucleotide (NAD), in its oxidized (NAD⁺) and reduced (NADH) forms. During glycolysis and the tricarboxylic acid cycle, NAD⁺ is reduced to NADH, while NADH is oxidized back to NAD⁺ in the mitochondrial respiratory chain and through lactic acid fermentation. Since NADH is autofluorescent and NAD⁺ is not, NADH fluorescence can serve as a proxy for the cellular redox state. In practice, NADH fluorescence is difficult to isolate from that of nicotinamide adenine dinucleotide phosphate (NADPH), a coenzyme involved in a wide range of anabolic processes. Nevertheless, shifts in the combined fluorescence signal of NADH and NADPH – referred to as NAD(P)H – are generally thought to be dominated by the influence of the NADH in the brain (Kasischke et al., 2004). Fluorescence lifetime imaging (FLIM) enables differentiation between free and protein-bound NAD(P)H molecules based on their distinct fluorescence lifetimes. Free NAD(P)H exhibits a shorter lifetime – around 0.4 ns – while bound NAD(P)H has a longer lifetime, ranging from 1 to 5 ns (Lakowicz et al., 1992; Kolenc and Quinn, 2019; Datta et al., 2020). By analyzing the contributions of these distinct lifetimes to the overall fluorescence decay, FLIM provides a quantitative estimate of the relative concentrations of free and bound NAD(P)H, robust to intensity distortion caused by light scattering. The free-to-bound NAD(P)H ratio measured by FLIM has been linked to changes in the balance between glycolysis and oxidative phosphorylation (OXPHOS), with a higher proportion of free NAD(P)H indicating an increased reliance on glycolysis relative to OXPHOS, and vice versa. This relationship has been observed in vitro and ex vivo in various contexts, including metabolic perturbations (Zheng et al., 2010; Stringari et al., 2012b; Drozdowicz-Tomsia et al., 2014; Liu et al., 2018; Yang et al., 2021), cell differentiation (Guo et al., 2008; Stringari et al., 2012a) and carcinogenesis (Bird et al., 2005; Skala et al., 2007; Rück et al., 2014), underscoring the utility of FLIM in assessing cellular metabolic states.

In the present study, we propose to demonstrate that NAD(P)H FLIM can capture the metabolic heterogeneities and shifts underpinning memory formation in the Drosophila brain. In contrast with hypothesis-driven approaches that rely on metabolite-specific sensors, NAD(P)H autofluorescence provides a global readout of cellular metabolism. Furthermore, autofluorescence imaging can be applied in non-genetic experimental models where genetically encoded indicators cannot be introduced. As a label-free and minimally invasive technique, NAD(P)H FLIM holds great promise for profiling neuronal metabolism in vivo. Despite this potential, its application to the brain has been so far limited. It has been used in vivo in the rodent brain to monitor strong and widespread metabolic responses to anoxia (Yaseen et al., 2013) and pharmacological interventions (Yaseen et al., 2017; Gómez et al., 2018)). However, its capacity to detect more subtle metabolic events associated with cognitive processes remains to be established.

Here, we employed FLIM of NAD(P)H to map the in vivo metabolic profiles of Drosophila neurons across multiple scales – from a broad view of the central brain to individual KC subtypes. This approach allowed us to visualize the cerebral metabolic heterogeneity at different levels and to identify a memory-related metabolic shift within a specific KC subtype. Our findings contribute to the understanding of the fine-scale energy regulation taking place in the MBs and establish NAD(P)H FLIM as a promising framework for future investigations into the energy dynamics underlying cognition.

Results

Metabolic landscape of the central brain highlights tissue-specific signatures

For this novel application of FLIM, we began by mapping the metabolic landscape across the entire Drosophila central brain, with a particular focus on the mushroom bodies (MBs). We conducted in vivo two-photon imaging of NAD(P)H FLIM over the entire central brain in flies expressing a DsRed cytosolic marker in MB neurons (MB-DsRed; Dataset 1 in Methods). Multi-exponential fitting and phasor analysis constitute the main approaches for NAD(P)H FLIM decay characterization (Digma et al., 2008; Stringari et al., 2011; Kolenc and Quinn, 2019; Schaefer et al., 2019; Datta et al., 2020). The former was chosen to automatically extract quantitative metabolic indicators. Given the novelty of the application, we opted for a simple and well-grounded approach to modeling NAD(P)H decays: the bi-exponential model (Kolenc and Quinn, 2019; Schaefer et al., 2019). The decay curves were spatially binned to enhance the signal-to-noise ratio (Figure 1a). A bi-exponential model was fitted to each decay using maximum-likelihood estimation (Sup. Figure 1), generating local estimates of the free NAD(P)H fraction – f_free – and the protein-bound NAD(P)H lifetime – τ_bound – across the central brain, for each fly. The weighted mean lifetime – τ_mean – was subsequently computed from these estimates.

The 3D spatial distributions of NAD(P)H properties could be visualized over the volume of the central brain (Figure 1b and Sup. Movies). This mapping revealed the distinct metabolic fingerprints of large anatomical structures. In particular, the neuropil could be well distinguished from the cortex region, which contains neuronal somata and cortex glia, as it displayed higher f_free and lower τ_bound values. Within the neuropil, layers with lower f_free and higher τ_bound values stood out, delineating different regions, as visible on the τ_mean map. Based on qualitative visual evaluation, these structures likely correspond to the layers formed by processes of neuropil- and tract-ensheathing glia, as well as tracheal processes (Pereanu et al., 2007; Kremer et al., 2017).

To quantitatively compare the different MB regions, we coarsely segmented the labeled MB image into defined anatomical structures using thresholding of the MB-specific fluorescent red marker and manually defined separation planes. Due to the limited resolution of the cytosolic marker image at this scale, the cortical region of the MB, primarily consisting of KC somata enveloped by cortex glia, could not be isolated from the adjacent calyx, which contains KC neurites and axonal projections from external neurons. The resulting segmentation included the peduncle, the medial lobes and the vertical lobe and the posterior part of the MB, consisting of the cortical region and calyx (Figure 1d).

The obtained segmentation enabled region-specific extraction of FLIM voxels, followed by the estimation of corresponding decay parameters. f_free and τ_bound distributions were predominantly unimodal (Sup. Figure 2), allowing for effective representation by the mean. No significant variation in the average f_free was observed across regions (Figure 1e, left). However, the posterior region showed notably higher τ_bound values than the rest of the MB (Figure 1e, right), consistent with what was observed at the whole brain level.

Overall, these results highlight regional differences in f_free and τ_bound among the neuropil, the cortex region and ensheathing glia in the central brain. At the MB level, a difference in τ_bound was still observed between the posterior part of the MB and purely synaptic regions. The interpretation of f_free as an indicator of the glycolysis-to-OXPHOS ratio has theoretical support (Schaefer et al., 2019) and has been validated in various contexts (see Introduction). The interpretation of τ_bound is more complex, as it is generally thought to represent a weighted average of fluorescence lifetimes of different enzyme-bound NAD(P)H species and to be influenced by changes in the composition of the proteins binding to NAD(P)H molecules (Kasischke et al., 2004; Yaseen et al., 2013, 2017; Evers et al., 2018; Sharick et al., 2018). Some studies also suggested that τ_bound variations may more sensitively reflect changes in the NADPH/NADH ratio than in NAD(P)H binding state, with τ_bound being positively correlated to the NADPH/NADH ratio (Blacker et al., 2014; Blacker and Duchen, 2016). Here, the difference between neuropil and cortical regions seems to be more strongly reflected by τ_bound than f_free, as confirmed by the quantitative analysis in the MBs, suggesting differences between these tissues in the composition of the proteins binding to NAD(P)H and/or in the NADPH/NADH ratio.

Figure 1

Somata and calyx regions have different metabolic profiles

The neuropil region of the MB – comprising the calyx, peduncle and lobes – contains KC neurites and is innervated by numerous additional neuronal types, including projection neurons, dopaminergic neurons, and MB output neurons as well as broadly-innervating DPM and APL neurons, which complicates the attribution of observed metabolic signals to particular cells. Isolating the NAD(P)H signal from the MB cortical region thus allows for the specific characterization of KC somata’s metabolic profile and its comparison with that of the adjacent synaptic region.

To this end, we acquired images with a more restrained scope, centered on the cortical area of the MBs. Each acquisition captured the KC-specific cytosolic marker MB-DsRed and NAD(P)H FLIM in a single hemisphere (Dataset 2 in Methods). We analyzed these images using mapping tools derived from other datasets: a MB template (Sup. Figure 3) and two segmentation masks, designed to isolate the KC somata and the calyx, respectively (Sup. Figure 4). The template image, designed to represent the average shape of the posterior MB region (Avants et al., 2010), was derived from a set of 165 MB-DsRed images (Dataset 3 in Methods) using a 3D volume registration procedure consisting of affine and non linear diffeomorphic transformations (Avants et al., 2008). A common template was established for left and right hemisphere MBs after confirming that their average shape differences were minimal (Sup. Figure 5). The masks were derived from MB-DsRed images and KC-specific nuclear marker images (Dataset 3 and Dataset 4 in Methods), using the same registration method. Once registered to each image, the masks enabled us to select the regions of interest in the FLIM image (Figure 2a). Decay parameters were then extracted as previously.

Registration was further used to project decay parameter images to the template space, enabling the generation of an average map of f_free and τ_bound in the cortical region and calyx (Figure 2b and Sup. Figure 6a) while minimizing the impact of inter-individual anatomical variations. This map was obtained by pooling left and mirrored right hemisphere images, after confirming that there was no hemisphere-specific spatial distribution of f_free (see Supplementary information and Sup. Figure 7). This average map highlighted that f_free was consistently higher in the calyx neuropil than in the cortical region of the MB.

To quantitatively substantiate this observation, the within-subject average was computed using both hemisphere images when available, or a single hemisphere image otherwise. The results confirmed that the calyx exhibits significantly higher values of f_free than the KC somata region, with means of 0.75 and 0.68, respectively (Figure 2b, bottom right). According to the literature, the difference in f_free likely reflects a higher reliance on glycolysis relative to OXPHOS in the calyx. Similar analysis of τ_bound revealed a significantly shorter bound-NAD(P)H lifetime in the calyx compared to the somata region (Sup. Figure 6b). This difference in τ_bound is consistent with the previously observed trends between the posterior region and the rest of the MB and, at a larger scale, between cortex and neuropil in general. Likewise, it could reflect differences in protein composition and/or in the NADPH/NADH ratio.

Measures of f_free and τ_bound were notably anticorrelated in the posterior part of the MB (Sup. Figure 8). The anticorrelation between these parameters can be qualitatively observed in several studies (Bird et al., 2005; Skala et al., 2007; Wetzker and Reinhardt, 2019; Yang et al., 2021). Multi-exponential fitting has been reported to potentially introduce artificial correlations between parameter estimates (Grinvald and Steinberg, 1974; Johnson, 2000). Through simulation, we confirmed that the relationship observed in our data was not an artifact of the fitting process (see Supplementary information and Sup. Figure 9), supporting the hypothesis of its biological origin.

Figure 2

Memory formation does not induce a significant metabolic shift over the MB somata and calyx regions

After imaging basal metabolism in the posterior MB region, we investigated whether memory formation induced detectable metabolic changes To this end, we acquired images of the MB posterior region in conditioned flies (refer to Dataset 5 in Methods). Flies were subjected to a classical olfactory aversive learning protocol (Tully and Quinn, 1985). In this paradigm, flies are exposed simultaneously to an odorant (odor A) and a train of electric footshocks, and then to another odorant without shocks (odor B; Figure 2c, top). Following one such training session, it is known that flies develop short-(STM) and middle-term memory (MTM) in the form of learned avoidance of odor A compared to odor B (Margulies et al., 2005). Control flies were subjected to unpaired conditioning, consisting of an exposure to the aversive stimulus and then to both odors, sequentially. No significant difference of the within-subject averages of f_free were observed in the somata or calyx regions between conditioned and control flies (Figure 2c, bottom). Similarly, no memory-associated shifts in τ_bound were detected in either region (Sup. Figure 6c).

In this experiment, indicators derived from NAD(P)H FLIM did not reveal a significant metabolic shift affecting globally the KC somata region or the calyx after a single training cycle. This result suggests the absence of an overall metabolic shift in these regions or its low amplitude relative to NAD(P)H FLIM sensitivity. Alternatively, specific effects may be diluted or averaged out, as the analyzed areas encompass diverse cell types that could exhibit distinct memory-related metabolic dynamics.

KC subtypes have different basal metabolic profiles

To further explore the metabolic heterogeneity of the MB, we aimed at characterizing the basal metabolism of the main subtypes of KCs (α/β, α’/β’ and γ). These subpopulations are known to have distinct functional roles. In particular, α/β neurons are specifically involved in middle- and long-term memories, whose formation has been shown to critically depend on energy regulation, while γ neurons are more specifically involved in short-term memory.

We initially attempted to separate subtype contributions within images through a mapping approach. To this end, we computed the average spatial distribution of each of these subtypes within the somata region of the template image (Sup. Figure 10a). While anatomical maps of the MB neuropil have shown that the neurites of different KC subtypes form distinct bundles (Aso et al., 2014; Li et al., 2020), the stereotypy of KC somata distribution remained largely unexplored. Our analysis revealed that these distributions were spatially non-uniform and reproducible from one group of flies to another (see Supplementary information and Sup. Figure 10bc). Comparisons of pairs of images from identical or different hemispheres revealed an absence of lateralization of the spatial distribution of the neuronal subtypes (see Supplementary information and Sup. Figure 11). However, the degree of overlap between the average subtype distributions was high (Sup. Figure 10de), suggesting important inter-individual variations and making effective separation through segmentation masks impossible.

We therefore performed recording in flies, expressing a nuclear marker in either α/β, α’/β’ or γ neurons, to separately measure the metabolic profile of each subpopulation while minimizing the contribution of the surrounding tissues. For each fly, standard two-photon imaging of this marker was performed alongside NAD(P)H FLIM over the somata region of the MB, one hemisphere at a time (Dataset 6 in Methods). Nuclear marker images were used to isolate the fluorescence decays around the somata of the marked neuronal subtype (Figure 3a). Within-subject averages of f_free were then computed from the masked FLIM images. These measurements revealed that the KC subtypes exhibit different basal metabolic profiles. Indeed, α/β, α’/β’ and γ somata displayed decreasing f_free values of 0.73, 0.71 and 0.70, respectively, with statistical testing showing a significant difference between α/β and γ subtypes (Figure 3b). Considering τ_bound, no statistically significant difference was observed between the three subtypes (Sup. Figure 12a).

Contrarily from previously observed regional differences, KC subtypes differed on f_free but not on τ_bound. This suggests, in line with the established interpretation of f_free, that the basal metabolism of α/β neurons rely more heavily on glycolysis relative to OXPHOS than that of γ neurons. We conducted an initial investigation to determine whether the lower f_free values observed in γ neurons could be attributed to a higher reliance on lactate-to-pyruvate conversion catalyzed by lactate dehydrogenase (Ldh) rather than pyruvate production by glycolysis, as Ldh was recently shown to be preferentially expressed in γ KCs as compared to other KC subtypes (Raun et al., 2023). However, we found no significant differences in f_free between Ldh-knockdown flies and controls (Dataset 7 in Methods and Sup. Figure 13).

This experiment provides initial insights into the metabolic heterogeneity of KC subtypes. Due to the limitations of the segmentation, it should be noted that our observations of the KC somata of a given subtype could partially include the surrounding glial processes.

Figure 3

Memory formation induces a metabolic shift in α/β KCs

After establishing the basal metabolic profiles of KC subtypes, we sought to measure the influence of memory formation on specific subtypes. We imaged the MB cortical region in conditioned flies expressing subtype-specific nuclear markers (Dataset 8 and Dataset 9 in Methods). We focused on α/β and γ neurons, which have been reported to support the formation of MTM and STM, respectively, following a single training cycle. Our analysis showed that flies exposed to paired training exhibited significantly lower f_free values compared to controls around α/β somata (Figure 3c). In contrast, the average f_free value in γ somata appeared unaffected by the type of conditioning. The effect of paired conditioning on τ_bound around α/β somata was noticeable but not statistically significant (Sup. Figure 12b).

These experiments demonstrated the capability of NAD(P)H FLIM to detect subtle physiological alterations of cellular metabolism induced by memory formation. Moreover, the results indicate that distinct metabolic dynamics unfold in α/β and γ neurons after aversive conditioning. More specifically, the decrease of f_free measured for α/β neurons suggests that, during memory formation, ATP is increasingly produced through OXPHOS relative to glycolysis.

Discussion

In this study, we sought to demonstrate the potential of NAD(P)H FLIM for detecting physiologically relevant spatiotemporal dynamics in brain energy metabolism. By applying this technique to the central brain of Drosophila melanogaster, we uncovered new insights into cellular metabolic heterogeneity. To establish this novel application, we employed bi-exponential fitting to extract two key parameters: the free fraction of NAD(P)H (f_free) and the lifetime of enzyme-bound NAD(P)H (τ_bound). Using genetically encoded anatomical markers combined with registration tools, these parameters were measured across different brain regions.

On a broad scale, our findings revealed the distinct metabolic profiles of ensheathing glial layers, cortical areas and neuropil regions, with τ_bound being particularly discriminative of the neuropil and cortex. At the MB level, differences between KC somata and calyx regions were observed in both τ_bound and f_free. At a finer scale, significant basal and memory-induced differences were detected in f_free among KC subtypes, whereas changes in τ_bound remained limited. These patterns suggest that tissue-level differences are primarily reflected by τ_bound, whereas subtype-specific variations among KCs are predominantly captured by f_free.

As FLIM-based metabolic measurements are influenced by numerous factors (Schaefer et al., 2019), the interpretation of these parameters in this novel context remains multifaceted. While f_free is often considered a reliable indicator of the glycolysis-to-OXPHOS balance, τ_bound is thought to reflect a broader range of factors, including the composition of proteins binding to NAD(P)H and the NADPH/NADH ratio. The tissue-level differences dominated by τ_bound, such as those observed between the posterior MB and other regions (Figure 1e), likely reflect broad metabolic distinctions that extend beyond energy-related processes. In contrast, the differences in f_free between KC subtypes may more reliably indicate distinct reliance on glycolysis versus OXPHOS. Instances where both f_free and τ_bound differ strongly – such as between the KC somata and calyx (Figure 2b, Sup. Figure 6a and Sup. Figure 8b) – might result from a combination of these factors. Measures of f_free and τ_bound in the posterior part of the MB and within KC subtypes were notably anticorrelated (Sup. Figure 8). We found, via simulation, that this anticorrelation was not an artifact of the fitting process, supporting the existence of a biological dependence between both parameters in our data (see Supplementary information and Sup. Figure 9).

At the level of KC subtypes, we found that α/β neurons exhibited higher f_free values than γ neurons, suggesting a higher glycolysis-to-OXPHOS ratio in α/β neurons. Ldh has been reported to be more expressed in γ neurons compared to the other KC subtypes (Raun et al., 2023), leading us to hypothesize that lactate-to-pyruvate conversion by Ldh reduces the glycolytic demand in this subtype. A knockdown of Ldh in γ neurons was conducted as a preliminary investigation but no significant impact on f_free was observed. Disentangling the contributions of KC somata and surrounding cortex glia to the FLIM signal, potentially through co-labeling approaches, could deepen our understanding of the basal metabolic heterogeneity among KC subtypes.

Notably, we observed a decrease in f_free in α/β somata – and surrounding glial processes – following a single aversive conditioning session, with no corresponding change in γ neurons. This reduction may reflect a transition towards OXPHOS in α/β somata during memory formation or maintenance. Although single-session memory training primarily induces the encoding of labile STM in γ neurons, parallel processes in α/β neurons lead to the formation of longer-persisting MTM (Blum et al., 2009; Séjourné et al., 2011; Bouzaiane et al., 2015). Our observation of a metabolic shift specific to α/β neurons aligns with recent evidence showing increased mitochondrial activity in KCs during MTM formation, supported by alanine transfer from glia (Rabah et al., 2023).

While interpretations of the observed NAD(P)H state variations in this new context remain open-ended, our findings lay a solid foundation for developing hypotheses about energy regulation in the MBs. This novel application of NAD(P)H FLIM can be extended in numerous ways, to probe the metabolic dynamics of memory formation across different conditions, resolutions and timescales. This approach also holds promise for investigating the role of metabolic dysregulation in age-related cognitive decline and neurodegenerative diseases (Saitoe et al., 2005; Camandola and Mattson, 2017). More broadly, it offers a powerful tool to advance fundamental neuroenergetic research, deepening our understanding of the interplay between energy metabolism and brain function.

Material and methods

Fly strains

Flies (Drosophila melanogaster) were raised on standard fly food (inactivated yeast 6% w/v; corn flour 6.66 % w/v; agar 0.9% w/v; methyl 4-hydroxybenzoate 22mM) on a 12-h light/dark cycle at 18°C with 60% humidity (unless mentioned otherwise). All experiments in this study were performed on 1-4-day-old adult female flies. Female flies were prefered for imaging experiments because of their larger size. The following split-GAL4 driver lines (Aso et al., 2014) were used to target specific populations of KCs: MB010B (all KCs; RRID:BDSC_68293), MB008B (α/β KCs; RRID:BDSC_68291), MB005B (α’/β’ KCs; RRID:BDSC_68306), MB009B (γ KCs; RRID:BDSC_68292). Reporter lines were UAS-Stinger (RRID:BDSC_84277), UAS-mCherry-NLS (RRID:BDSC_38424; Delestro et al., 2020), and UAS-Ldh^RNAi (HMS00039; RRID:BDSC_33640), as well as MB-DsRed (Riemensperger et al., 2005).

Olfactory conditioning

Olfactory conditioning was conducted similarly to other studies from the same research group (Pascual and Préat, 2001; Isabel et al., 2004; Plaçais et al., 2012; Rabah et al., 2023). Experimental flies were transferred to fresh bottles containing standard medium on the day before conditioning. Groups of ∼30 flies were subjected to a single cycle olfactory conditioning protocol, using either paired or unpaired (control) stimuli in custom-made barrel-type machines (Pascual and Préat, 2001). Throughout the conditioning protocol, each barrel was attached to a constant air flow at 2 l/min. For a single cycle of associative training, flies were first exposed to an odorant A for 1 min while twelve 60-Volt 1.25-second square-wave pulses were delivered at 5-second intervals; flies were then exposed 45 s later to a second odorant B without shocks for 1 min (Figure 2c). The odorants 3-octanol and 4-methylcyclohexanol, diluted in paraffin oil at a final concentration of 2.79·10^-1 g.l^-1, were used as odorants A and B, respectively. During unpaired conditionings, the shock and odor stimuli were delivered sequentially, with the first odor being released 2 min after the end of the shocks.

Image acquisition

Data were collected indiscriminately from 30 min to 1.5 h after memory training. A single fly was picked and prepared for imaging as previously described (Plaçais et al., 2017). Briefly, the fly was glued on a plastic slide using a biocompatible dental glue (3 M ESPE Protemp) without any prior anesthesia. An alignment wire was used to maintain the Drosophila head in a correct position. The head capsule was opened and the brain was exposed by gently removing the superior tracheae. The head capsule was bathed in artificial hemolymph solution for the duration of the preparation. The composition of this solution was NaCl 130 mM (Sigma cat. no. S9625), KCl 5 mM (Sigma cat. no. P3911), MgCl2 2 mM (Sigma cat. no. M9272), CaCl2 2 mM (Sigma cat. no. C3881), D-trehalose 5 mM (Sigma T cat. no. 9531), sucrose 30 mM (Sigma cat. no. S9378) and HEPES hemisodium salt 5 mM (Sigma cat. no. H7637). At the end of surgery, any remaining solution was absorbed and a fresh 90-μl droplet of this solution was applied on top of the brain.

Two-photon fluorescence lifetime measurements were performed using a Leica Microsystems SP8 DIVE-FALCON microscope equipped with a HC-IRAPO 25x, 1.0 NA water-immersion objective coupled to a dual output InSight infrared excitation laser (Spectra Physics). Two-photon excitation of NAD(P)H autofluorescence was achieved using the tunable laser beam at 740 nm. Fluorescence was measured by time-correlated single-photon counting in the spectral range 420-500 nm. Fluorescent anatomical markers were imaged simultaneously. Ds-Red and mCherry were excited at 740 nm and their emission intensity was recorded in the spectral range 600-640 nm. GFP was excited at 910 nm and its emission intensity was recorded in the spectral range 500-550 nm. All data acquisition and analysis was performed with the Leica LAS-X (v3.5.7) software. 3D z-stacks were acquired either on the whole brain, or on the cortex region encompassing all Kenyon cells. The number of planes in the z-dimension was adjusted to include all visible KCs somata. The excitation laser intensity was adjusted with LAS-X controls to 14% corresponding to ∼9 mW at the objective output. On each x-y plane, 30 sequential scans were averaged, yielding on average 150-300 photons per pixel.

Datasets

This section describes the different image datasets used throughout the present study.

Dataset 1 consists of 5 images from 5 flies expressing a KC-specific cytosolic marker (MB-DsRed). This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the entire central brain, with an average resolution of 0.74 µm for the X-Y axes and 2.00 µm for the Z axis.

Dataset 2 consists of 31 images from 17 flies expressing a KC-specific cytosolic marker (MB-DsRed).. This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.40 µm for the X-Y axes and 2.00 µm for the Z axis. For additional details, see Sup. Table 1.

Dataset 3 consists of 165 images from 102 flies expressing a KC-specific cytosolic marker (MB-DsRed) and a subtype-specific nuclear marker (MB008B>mCherry::NLS; MB009B>mCherry::NLS; MB005B>mCherry::NLS). Both markers were simultaneously imaged in separate channels. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.35 µm for the X-Y axes and 0.60 µm for the Z axis. For additional details, see Sup. Table 2.

Dataset 4 consists of 32 images from 19 flies expressing a KC-specific nuclear (MB010B>Stinger) and cytosolic (MB-DsRed) markers. Both markers were simultaneously imaged in separate channels. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.38 µm for the X-Y axes and 0.60 µm for the Z axis. For additional details, see Sup. Table 3.

Dataset 5 consists of 31 images from 31 conditioned flies expressing a KC-specific cytosolic marker (MB-DsRed).. This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.39 µm for the X-Y axes and 2.00 µm for the Z axis. For additional details, see Sup. Table 4.

Dataset 6 consists of 79 images from 43 flies expressing a subtype-specific nuclear marker (MB008B>mCherry::NLS; MB009B>mCherry::NLS; MB005B>mCherry::NLS). This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.38 µm for the X-Y axes and 4.00 µm for the Z axis. For additional details, see Sup. Table 5.

Dataset 7 consists of 24 images from 24 flies, with and without Ldh-knockdown mutations (MB009B>Ldh RNAi), expressing a γ-specific nuclear marker (MB009B>mCherry::NLS). This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.37 µm for the X-Y axes and 2.00 µm for the Z axis. For additional details, see Sup. Table 6.

Dataset 8 consists of 40 images from 38 flies expressing an α/β-specific nuclear marker (MB008B>mCherry::NLSr). This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.38 µm for the X-Y axes and 2.00 µm for the Z axis. For additional details, see Sup. Table 7.

Dataset 9 consists of 23 images from 23 flies expressing a γ-specific nuclear marker (MB009B>mCherry::NLS). This marker was imaged simultaneously with NAD(P)H FLIM. Each 3D image captured the somata and calyx region of a single hemisphere, with an average resolution of 0.39 µm for the X-Y axes and 2.00 µm for the Z axis. For additional details, see Sup. Table 8.

Segmentation of fluorescence intensity images for FLIM masking

For Dataset 1, each image was initially divided into two halves along a manually defined median plane, enabling separate processing of both hemispheres. The somata and calyx region was roughly isolated from each image using manually defined planes orthogonal to the Y axis. A binary mask of this region was then obtained through Li thresholding (Li and Tam, 1998) after Gaussian smoothing (with a standard deviation σ = 4 µm). In each single-hemisphere image, the entire MB was delineated by applying Li thresholding after Gaussian smoothing (σ = 2 µm). The peduncle and lobes were first roughly isolated using manually defined planes, orthogonal to the main axes. Binary masks for these regions were then created using smoothing (σ = 2 µm) and Otsu thresholding (Otsu, 1979). Finally, all single-hemisphere masks were combined to reconstruct complete brain masks.

Images of Dataset 2 and Dataset 5 were segmented by applying Li thresholding after Gaussian smoothing (σ = 2 µm) and keeping only the largest contiguous volume.

Images of Dataset 6, Dataset 7, Dataset 8 and Dataset 9 were segmented in two steps. First, Li thresholding was applied after Gaussian smoothing (σ = 0.5 µm) and only the volumes superior to 1 µm³ were kept. Then, Li thresholding was applied again within the resulting non-zero voxels.

FLIM image processing

We observed that voxels with exceptionally high photon counts were predominantly associated with tracheas or with objects outside the region of interest due to the limitations of the segmentations. To reduce the influence of these artifacts while preserving the neuronal signal, we excluded voxels with photon counts exceeding Q₃ + 3 × IQR, where Q₃ represents the upper quartile of the count distribution and IQR denotes the interquartile range.

To ensure that the total photon counts of the decay histograms were sufficient for reliable curve fitting, the signals of neighboring voxels were combined. This process, known as spatial binning, enhances the signal-to-noise ratio of the decay signal at the cost of reducing spatial resolution. Decay histograms were summed within a disk centered on each voxel. This disk was defined as a 5 x 5 x 1 voxel volume, with the four corner voxels excluded. To avoid unreliable modeling of the decays, voxels for which the total photon count did not exceed 500 after binning were excluded from further analyses.

Theoretical fluorescence decays were modeled using the following bi-exponential function D:

where f_free is the free NAD(P)H fraction, τ_freeis the free NAD(P)H lifetime and τ_bound is the protein-bound NAD(P)H lifetime. The measured decays were fitted with the following function f, which incorporates the theoretical decay and the effects of the measurement process:

where IRF is the microscope impulse response function, Δt is the temporal shift between the IRF and the measured signal and T is the laser period. The modeled signal at time t takes into account the fluorescence caused by the last laser pulse (term for p = 0) as well as the previous one (term for p = 1), as proposed by (Klemm et al., 2015). (Klemm et al., 2015).

In practice, f was evaluated for N_t = 129 time steps t_i corresponding to the time bins of the FLIM data and normalized:

The final simulated decay sequence s was defined as:

where b is a constant background noise.

The IRF was obtained from each image using the FALCON module of the Leica LAS X software managing FLIM acquisitions (v3.5.7). T was set to 12.5 ns, corresponding to the 80 MHz repetition rate of the microscope laser. Due to its limited physiological interpretability and its narrow reported range of variation (Guo et al., 2008), τ_freewas fixed at 0.4 ns, in line with values reported in the literature for aqueous solutions (Scott et al., 1970; Visser and van Hoek, 1981), cell cultures (Wu et al., 2006; Rück et al., 2014; Sharick et al., 2018), tissue slices (Vishwasrao et al., 2005), in vivo experiments (Yaseen et al., 2013), as well as in our own findings (Sup. Figure 14).

Parameters f_free, τ_bound, Δt and b were estimated via maximum-likelihood estimation (MLE), using the 2I* quality parameter (Maus et al., 2001):

where m_i represents the measured photon count in bin i and C represents the total photon count of the measured decay (i.e the sum of m_i values). MLE was preferred over the least squares method for its higher reported accuracy (Maus et al., 2001; Laurence and Chromy, 2010; Chen et al., 2022).

This fitting process was implemented using the Nelder-Mead algorithm (Gao and Han, 2012). Each curve fitting was repeated 5 times, using initialization parameters uniformly sampled from empirically defined ranges. The parameter estimation having returned the lowest 2I* value was selected. Voxels with final 2I* values exceeding the empirically defined threshold of 200 were systematically excluded from further analysis (Sup. Figure 1).

To accelerate the estimation for all voxels of an image, a two-step process was carried out. An initial fitting was computed on a downsampled version of the image, which had been decimated by a factor of 3 in the X and Y directions. The time shift Δt was then fixed to the median of the obtained values, enabling the application of a streamlined three-parameter model to the complete image.

In the end, we considered that the metabolic profile extracted from a given decay curve is constituted of the couple (f_free, τ_bound). The mean lifetime was optionally computed as the sum of the free and bound NAD(P)H lifetimes weighted by their fractional contributions:

Analysis of the spatial distributions revealed that voxels with f_freevalues below 0.25 consistently fell outside the regions of interest as a result of imperfect segmentation. Consequently, these voxels were systematically excluded from further analysis.

Building a template image of the somata and calyx region

A template image of the somata region of the MB was generated from the 165 KC-specific cytosolic marker channel of the images (Dataset 3). The preprocessing steps for these images are outlined in Sup. Figure 3a. Images were resampled to an isotropic voxel size of 0.5 µm in all directions. To eliminate background noise and extraneous bright objects, we applied Li thresholding to a smoothed version of each image, retaining only the brightest contiguous object. The images were cropped and padded to make their size uniform and to center them based on their center of mass. The average intensity histogram of the dataset was computed by averaging the intensity histograms of all images. The intensity of all images were normalized by adjusting their intensity histogram to match the dataset’s average histogram (Nyúl et al., 2000).

The template was subsequently constructed from the preprocessed images using the symmetric group-wise normalization algorithm (Avants et al., 2010). The registration method used in the procedure was defined as a sequence of rigid, affine and diffeomorphic transformations applied to the 3D volumetric images. Symmetric normalization (Avants et al., 2008) was used as the diffeomorphic transformation, with empirically defined parameters of 0.2 for the gradient step, 3 for the update field smoothing and 1 for the displacement field smoothing. The method initializes the template as the mean image and then iteratively updates it to minimize the average displacement required to register the dataset to the template (Sup. Figure 3b). The template image associated with the lowest average displacement after 8 iterations was selected. The template was finally refined by applying thresholding and centering. The final template image is shown in Sup. Figure 3cd.

Building masks for the somata and calyx region

A mask isolating the somata region was derived from the images of Dataset 4, obtained by simultaneously capturing KC-specific nuclear and cytosolic markers on two separate channels. This set of images was registered to the previously established MB template (Sup. Figure 4) using the cytosolic marker channel. The nuclear marker channels were warped using the obtained transformations and then averaged. The somata region mask was obtained by binarizing the result. The binarization used the highest threshold obtained through the application of the 3-class multi-Otsu algorithm (Liao et al., 2001) after Gaussian smoothing (σ = 3 µm).

A mask isolating the calyx region was derived by binarizing the MB template image. The binarization used the highest threshold obtained through the application of the 3-class multi-Otsu algorithm after Gaussian smoothing (σ = 5 µm). The final calyx mask was refined by excluding the voxels included in the somata mask.

Mapping f_free in the somata and calyx region

Single-hemisphere maps of f_free were obtained from Dataset 2. Right-side images were mirrored to allow their superposition with left hemisphere images. The MB template was registered to each image. This transformation was applied to the chosen mask (somata or calyx) using nearest-neighbor interpolation. The warped mask was applied to select the region of interest in the FLIM image. Spatial binning and curve fitting was applied to the masked FLIM image, resulting in maps of the decay parameters over the region of interest.

To map decay parameters in the template space, the KC-specific marker channel (MB-DsRed) of the image was registered to the template. This deformation was applied to the maps of decay parameters using nearest-neighbor interpolation. The average f_free map presented in Figure 2b was obtained by averaging all registered parameters maps.

Statistics

For all boxplots, whiskers show the minimum and maximum values, the box shows the first and third quartiles and the horizontal line displays the median. All tested distributions are described by their sample size (n = sample size) and their mean and standard deviation (µ = mean ± standard deviation). All Student’s t-tests are two-tailed independent samples t-tests, unless stated otherwise. In the case of multiple comparisons, displayed and reported p-values are the Bonferroni-adjusted p-values, as specified in the corresponding captions.

Data availability

Data, including FLIM, marker imaging and IRFs, is available in the BioStudies database (https://www.ebi.ac.uk/biostudies) under accession number S-BIAD1528 (DOI: 10.6019/S-BIAD1528).

Code availability

The code used to process marker and FLIM images is available as a GitHub repository (https://github.com/biocompibens/MB-FLIM).

Acknowledgements

This research has received support from the Region Ile de France through the DIM Elicit program, from the Institut Convergences QLife (ANR-17-QLIFE), from the LabEx MemoLife (ANR-10-LABX-54 MEMOLIFE and ANR-10-IDEX-0001-02 PSL* Université Paris), from the European Research Council (ERC-AdG-741550), from the Human Frontier Science Program (RGY0078/2017 ChroMet) and from the Agence Nationale de la Recherche (ANR) under contracts ANR-10-INBS-04 France BioImaging and ANR-11-EQPX-0029 Morphoscope2.

Supplementary information

Supplementary results

The template image minimizes displacement upon registration

We confirmed that registering the images to the template image resulted in less displacement compared to registration to either the average image or to randomly selected images (Sup. Figure 3ef). This finding supports the use of the template image as a reference space to minimize deformation during registration.

Minimal differences in hemisphere-specific MB shapes

We compared templates generated from left, right or both hemispheres images of the somata region of the MBs (Sup. Figure 5a). Differences between left- and right-specific template images were consistently higher than differences obtained using mixed hemispheres (Sup. Figure 5b). As the subset size increased, this effect became stronger. For subsets of 30 images, the average differences were 11.7% and 14.4% for inter-hemisphere and mixed-hemisphere comparisons, respectively. This result suggests the existence of a slight difference in shape between left and right MBs. We concluded upon visual inspection that this putative difference was too subtle to justify the use of separate templates (Sup. Figure 5cd).

The spatial distributions of KC subtype somata are not hemisphere-specific

We compared the spatial distributions of the somata of KC subtypes in images of Dataset 3. Using pairs of MB images coming from different subjects, the distances between left and right hemispheres were not significantly different from the intra-hemisphere differences (Sup. Figure 11). This result indicated that the subtype distributions are not hemisphere specific. This analysis also showed that distances between left and right α/β distributions were lower when they came from the same subject.

Somata of KC subtypes have specific spatial distributions

We tested whether the average spatial distribution of the somata of each KC subtype was specific, using images of Dataset 3 featuring a KC-specific marker channel and a subtype-specific marker channel. Image subsets of increasing size were used to compare the average spatial distributions of the different subtypes (Sup. Figure 10b). Average distributions from images of the same subtype consistently showed greater similarity than to those from different subtypes (Sup. Figure 10c). When using 20 images per spatial distribution, the average differences were 8.1%, 13.5% and 13.1% for intra-α/β, intra-α’/β’ and intra-γ comparisons, respectively, and 19.4% for inter-type comparisons. This result demonstrates that the spatial distribution of KC somata is subtype-specific.

The average spatial distributions of KC subtypes are highly overlapping

We have shown that the average spatial distributions of the somata of KC subtypes were different. These distributions can be visualized through the subtype proportion map (Sup. Figure 10d). As can be observed from the color scale of this map, the proportions of each subtype vary within specific ranges: 20-70%, 10-40% and 10-60% for α/β, α’/β’ and γ, respectively. When looking at the distributions of the proportions, we observe that in most voxels, the proportion of a given subtype corresponds to the overall fraction of neurons belonging to this subtype (Sup. Figure 10e). This shows that although the distributions are non-uniform, the high inter-subject variability prevents the delineation of areas that consistently isolate specific neuronal subtypes across individual flies.

The spatial distribution of f_free in the MB cortical region is not hemisphere-specific

We compared the spatial distributions of f_free in images of Dataset 2. Using pairs of MB images coming from different subjects, the distances between left and right hemispheres were not significantly different from the intra-hemisphere differences (Sup. Figure 7). This result indicated that the subtype distributions are not hemisphere specific.

The anticorrelation between f_free and τ_bound is not artificially introduced by the fitting process

It has been long known that estimations of parameters of multi-exponential models tend to be correlated, due to the fact that multiple parameters can affect the simulated curve in a similar way (Grinvald and Steinberg, 1974; Johnson, 2000). In the case of NAD(P)H FLIM, increasing f_free or reducing τ_bound affects the overall fluorescence decay similarly. We hypothesized that this artifact of the fitting process could partially explain why τ_bound and f_free values consistently exhibit anticorrelated patterns in multiple previous studies (Bird et al., 2005; Skala et al., 2007; Wetzker and Reinhardt, 2019; Yang et al., 2021) and in our data (Sup. Figure 8). To verify this hypothesis, we simulated 5,000 decay curves using our bi-exponential model, with τ_free set to 0.4 ns. b was fixed to 1.5·10^-3 to match the overall average of Dataset 2. τ_bound was uniformly sampled between 3 and 4.5 ns and f_free between 0.6 and 0.8, similarly to the ranges of within-subject averages obtained from Dataset 2. All simulations used the same arbitrarily selected IRF from Dataset 2 and Δt was set to 0. The decay curves were scaled to simulate a photon count of 6000, matching the overall average of Dataset 2. Finally, Poisson-distributed noise was applied. All simulated decay curves were fitted using the same bi-exponential model as in other analyses, with τ_bound and f_free as the only parameters. We observed a very low correlation of 0.06 between the estimated values of τ_bound and f_free (Sup. Figure 9a), implying that the anticorrelations observed in our study were unlikely introduced by the fitting process.

As a supplementary verification, we simulated 5,000 decay curves from correlated τ_bound and f_free couples. Correlated random variables C₁ and C₂ were obtained from independent random uniform variables X₁ and X₂using the Cholesky decomposition:

with the Pearson correlation ρ set to -0.86 to match the value obtained from Dataset 2. C₁and C₂ were then scaled and offset to vary in the same ranges as τ_bound and f_free in the previous simulation. In this simulation, the estimated parameters showed an important negative correlation of -0.64 (Sup. Figure 9b), confirming the capacity of the fitting process to retrieve dependencies between τ_bound and f_free.

Supplementary methods

Image distance metrics

To compare pairs of f_free spatial distributions registered in the template space, we used a metric termed the “f_free map distance”. This metric was defined as the average of absolute voxel-wise differences across all valid voxels in a pair of f_free maps. For the f_free maps of images A and B, noted F_A and F_B, this distance can be formulated as:

where N_mask is the number of voxels in the mask.

For fluorescence intensity images, we employed a metric termed the “normalized image distance”. The compared images were first normalized so that their voxel intensities summed to one. The normalized distance was then calculated as the sum of absolute voxel-wise differences between the normalized images, divided by 2. It can be understood as the fraction of the image intensity that should be displaced to match one image to the other. For images A and B, noted I_A and I_B, and their summed intensity S_A and S_B, this distance can be formulated as:

where N_image is the number of voxels in the images.

Building a map of subtype proportions

A map, estimating the average proportion of each KC subtype over the MB template was derived from Dataset 3. The images of this set are composed of 2 channels, one labeling all KCs with a cytosolic marker and the other labeling KCs from one KC subtype with a nuclear marker. Right-hemisphere images were mirrored. All images were registered to the template based on the KC-specific marker channel. We applied a previously established mask (Sup. Figure 4) to isolate the signal of the somata region in the subtype-specific channel. The masked images were normalized to make the total intensity of each image equal to the number of neurons of the associated subtypes (1002 α/β, 370 α’/β’ and 671 γ neurons on average according to Aso et al. (Aso et al., 2009). The intensity of the normalized images can be seen as an approximation of the neuronal density, expressed in neuron/voxel. The three subtype density maps are then summed to produce a KC density map. Finally each subtype density map is divided voxel-wise by the KC density map to produce a subtype proportion map.

Comparing hemisphere-specific MB shapes

This process was carried out using images of Dataset 3. All images were preprocessed following the same steps as for template building (Sup. Figure 3a). Our testing process involved repeated distance measurements (Sup. Figure 5a). The complete set of preprocessed images was first randomly divided into two non-overlapping subsets, arranged in random order. In the mixed-hemisphere case, the half sets were created from left and right images indiscriminately. In the inter-hemisphere case, the first and second subset exclusively contained left and right MB images, respectively. A template was constructed using the first N images from each subset. The resulting average spatial distributions were compared by computing the total absolute difference between them. This procedure was repeated multiple times, with 5 random splits performed for each case.

Testing the hemisphere-specificity of the spatial distributions of KC subtype somata

This process was carried out using images of Dataset 3. All images were registered and masked following the same steps as for the establishment of the subtype proportion map (Sup. Figure 10a). The normalized image distance was computed between all pairs of images. All distances were labeled as either intra-subject, inter-subject intra-hemisphere or inter-subject inter-hemisphere.

Testing the specificity of the spatial distributions of KC subtype somata

This process was carried out using images of Dataset 3. All images were registered, masked and normalized following the same steps as for the establishment of the subtype proportion map (Sup. Figure 10a). Our testing process involved repeated distance measurements (Sup. Figure 10b). For each measurement, the dataset was split into two non-overlapping subsets of images. Each set was composed of images displaying a unique subtype-specific marker. Intra-type comparisons were computed using two sets of the same subtype while extra-type comparisons used sets of different subtypes. In both cases, the first N images of each set were averaged. Finally, the total absolute difference was computed between the two resulting average spatial distributions. This scheme was repeated multiple times for different random splits and values of N. For each neuronal subtypes, 20 splits were used.

Testing the hemisphere-specificity of the spatial distribution of f_free in the somata region

This process was carried out using images of Dataset 2. Right-hemisphere images were mirrored. Maps of ffree in the somata region were obtained by registering the corresponding mask to each image. The resulting ffree maps were registered to the template using the KC-specific marker channel. The total absolute difference was computed between all pairs of registered f_free maps. All distances were labeled as intra-subject, inter-subject intra-hemisphere or inter-subject inter-hemisphere.

Supplementary figures

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Supplementary Figure 8.

Supplementary Figure 9.

Supplementary Figure 10.

Supplementary Figure 11.

Supplementary Figure 12.

Supplementary Figure 13.

Supplementary Figure 14.

Supplementary tables

Supplementary table 1

Supplementary table 2

Supplementary table 3

Supplementary table 4

Supplementary table 5

Supplementary table 6

Supplementary table 7

Supplementary table 8