Computational exploration of molecular receptive fields in the olfactory bulb reveals a glomerulus-centric chemical map

Progress in olfactory research is currently hampered by incomplete knowledge about chemical receptive ranges of primary receptors. Moreover, the chemical logic underlying the arrangement of computational units in the olfactory bulb has still not been resolved. We undertook a large-scale approach at characterising molecular receptive ranges (MRRs) of glomeruli in the dorsal olfactory bulb (dOB) innervated by the MOR18-2 olfactory receptor, also known as Olfr78, with human ortholog OR51E2. Guided by an iterative approach that combined biological screening and machine learning, we selected 214 odorants to characterise the response of MOR18-2 and its neighbouring glomeruli. We found that a combination of conventional physico-chemical and vibrational molecular descriptors performed best in predicting glomerular responses using nonlinear Support-Vector Regression. We also discovered several previously unknown odorants activating MOR18-2 glomeruli, and obtained detailed MRRs of MOR18-2 glomeruli and their neighbours. Our results confirm earlier findings that demonstrated tunotopy, that is, glomeruli with similar tuning curves tend to be located in spatial proximity in the dOB. In addition, our results indicate chemotopy, that is, a preference for glomeruli with similar physico-chemical MRR descriptions being located in spatial proximity. Together, these findings suggest the existence of a partial chemical map underlying glomerular arrangement in the dOB. Our methodology that combines machine learning and physiological measurements lights the way towards future high-throughput studies to deorphanise and characterise structure-activity relationships in olfaction.


Introduction
For most sensory systems, the receptive fields of primary receptor neurons are fairly well identified, for example the spectral tuning of photoreceptors, or frequency tuning of cochlear hair cells. Moreover, topographical arrangements like retinotopy in vision, tonotopy in audition, and somatotopy in somatosensation have been described [1][2][3] , which minimize wiring length and likely improve the efficiency of computation 4-6 . In the olfactory system, the organisation of the main olfactory epithelium (MOE) in zones that define preferred locations of olfactory sensory neuron (OSN) is an instance of such topography [7][8][9] . Some zonal co-localisation of OSNs with similar molecular receptive ranges (MRRs) has been observed [10][11][12] , but well-defined zonal adherence can only be observed for a subset of receptors. It has been suggested that zonal OSNs expression patterns propagate to the olfactory bulb (OB) 13 .
Conflicting interpretations of experimental results exist however as to whether glomerular positions in the OB reflect a chemical topology, in the sense that spatial proximity of glomeruli indicates proximity of their preferred ligands in chemical space. Mori et al. described a patchy map of clusters of neighbouring glomeruli with similar MRRs, located in stereotypic positions, as summarised in a review 14 along with other studies that suggest a "systematic, gradual, and multidimensional" representation of molecular features on the OB's surface. Similarly, Johnson and Leon reported that odorants sharing structural chemical features such as functional groups, aromatic rings and solubility evoked clustered responses 15 , that shift in a systematic way according to carbon chain length, confirming earlier observations 16 .
Contrarily, other reports dispute that chemical similarity is reflected in the bulbar location of glomeruli. Soucy et al. found that neighbouring glomeruli tend to respond to overlapping sets of ligands, but dismiss the idea of a chemotopic ordering, claiming that the intrinsic variability of glomerulus placement is greater than the accuracy of the chemical map 17 . They also shine a critical light on Mori's hypothesis on chemotopic modules as they observe glomeruli in a region apparently tuned to aldehydes that do not respond to that substance class at all. Generally, they could not confirm any chemotopic arrangement beyond a few glomerular radii.
Ma and colleagues 18 reported that no chemotopic organisation could be observed in a large dataset of glomerular responses in the dOB using about 60 odorants. However, their analysis mainly rests on functional group comparisons and largely ignores other molecular features, therefore their assessment does not reflect all aspects of chemical similarity. Their report was opposed by Yablonka and coworkers 19 , who pointed out flaws in their methodology regarding the normalisation of descriptors, and proposed an analysis that indeed suggests a relation between spatial distance of glomeruli and chemical distance of their preferred ligands.
These seemingly conflicting reports reveal the general problem in assessing chemotopy, namely that no generally valid, rigorous definition of chemical similarity exists. For example, similarity of odorants can be described by physico-chemical properties 20,21 . Likewise, molecular similarity can be assessed by vibrational spectra, which have been used to successfully predict receptor affinity, as originally proposed in drug discovery 22 , and more recently confirmed in insect olfaction 23 .
An important distinction to make when considering the seemingly contradictory evidence is whether actual chemoptopy has been assessed by calculating the overlap of MRRs in chemical space, or rather tunotopy, that is, merely computing the amount of overlap in the ligand spectra 18 . Tunotopy, while straightforward to compute empirically, lacks an immediate connection to physico-chemical space, thus preventing further insights into the chemical map on the OB surface. This chemical map requires an abstraction away from discrete compounds, and towards physico-chemical properties that define the map "coordinates". Due to the high dimensionality of chemical space, a reliable assessment of chemotopy thus critically depends on the number of odorants and glomerular responses taken into account.
MRR data are still sparse, and the sheer number of OR genes that are still to be de-orphaned and characterised has so far impaired a comprehensive view on chemotopy in the OB. Chemical tuning of a glomerulus and its neighbourhood is rarely assessed using more than a hundred odorants, which provides a tiny window into the vast and extremely high-dimensional olfactory chemical space. Detailed MRRs are only available for a small amount of receptors [24][25][26][27] , partial MRRs merely for a few more 20,[28][29][30][31][32][33] . The scarcity of data limits the explanatory power that can be expected from any approach that pursues the investigation of chemotopy organisation in the OB that relies on comparisons of a large number of glomerular MRRs.
In this study we address the scarcity of data by automated biological in-vivo screening to deorphanize glomeruli, characterise their MRRs, and we leveraging machine learning to assess tunotopic and chemotopic relations between neighboring glomeruli. We started by mapping the MRR of a single glomerulus as completely as possible, and subsequently explored the chemotopic embedding in its neighbourhood. This contrasts earlier attempts that third to investigate the topographic relationship of many glomeruli with fragmentary MRRs.
We initially focused on the class-1 receptor MOR18-2 (a.k.a. MOL2.3, or OLFR78) using mice with GFP labelled olfactory sensory neurons (OSNs), then explored MRRs of adjacent glomeruli, finally extending the analysis on the whole area covered by imaging, assessing chemotopy on several spatial scales. We used machine learning to obtain a physico-chemical activation model of MOR18-2, thus defining its MRR by means of physico-chemical properties.
This allowed us to map the location of a glomerulus in chemical space to its location on the dOB, getting further insight into chemotopy.

Methods
Mice were anaesthetized using urethane (1.5 g/kg i.p.). Anaesthetic was supplemented throughout the experiments and the body temperature was kept between 36.5C and 37.5C using a heating pad and a rectal probe. For imaging a craniotomy over one or both olfactory bulb was cut. The dura mater was removed and the imaging chamber was filled with agar (1.5 %) and covered with a glass cover slip. The prepared skull was fixated with cement to a metal plate under the microscope. All animal care and experimental procedures were approved by the regional authorities and carried out in accordance with the animal ethics guidelines of the Max Planck Society.
Odorants were filled into glass vials under Argon atmosphere in order to avoid any contamination. A two-arm robot (Combipal, CTC-Analytics, Zwingen, Switzerland) using the Software Chronos (Axel Semrau, Sprockhoevel, Germany) was used to sample 2.5 ml of the odour headspace and inject it into a constant carrier flow of filtered and humidified air (2 l/min) towards the mouse's nose. Odour molecules reached the nose 2.5 s ± 0.3 s after recording onset as measured by a photoionization detector (Aurora Scientific, Canada). We measured the response to different stimuli sets of 4 to 53 odours in 78 MOL2.3-IGITL mice 34 , 3 heterozygote OMP-SpH mice 35 and 15 heterozygote MOL2.3-IGITL OMP-Sph cross-breeds. Each stimulus set was at least measured twice in each mouse with odours presented in a pseudo-randomized sequence. We used mice of both sexes, aged 5-23 weeks.
After each odour presentation the syringe used for odour transfer was flushed with nitrogen for 72 s to minimize contamination. The sequence of odour presentation was pseudo-randomised to further quell any systematic influence of potential contamination. In addition, we recorded periods was presented at the mean point between two measurements. Those 'blank' periods were later used for bleaching correction in SpH-imaging (see below).

Intrinsic Optical Signal (IOS) Imaging
Odour responses were recorded in the dorsal olfactory bulb for 12s at 5Hz using a macroscope (Pentax zoom lens 12-48 mm, f = 1:1.0 and Nikkor lens 135 mm, f = 1:2.0) and an Orca-R2 camera (Hamamatsu, Japan; 1024 x 1344 px) under illumination with red light (690 nm) with a frame rate of 5 Hz (60 frames in 12 seconds). Uni-lateral OB recordings were performed at a focal length of 24 mm (field of view 1.63 mm x 1.24 mm), whereas bilateral recordings were performed at a focal length of 48 mm (field of view 3.26 mm x 2.48 mm). Before and after each presentation of the entire stimulus set, the pattern of blood vessels was recorded using green illumination (546 nm, 'green image') and controlled for shifts to exclude movement artefacts. Furthermore, the MOR18-2 glomeruli were located using blue illumination (475 nm, 'GPF image') and an emission filter at 535 nm.
To increase signal-to-noise ratio and reduce computational load we binned the raw data with an 8 x 8 px spatial and a 12 frame temporal window. Then the odour induced activation was calculated as the relative decrease of reflectance ∆ ⁄ = −( − ' )⁄ . ' was the mean reflectance during the first 2 s after recording onset, before the odours reached the nose (which happened at 2.5 s ± 0.3 s, see above). Furthermore the data was spatially bandpassfiltered with two Gaussian filters ( *+, = 10 px, 0120 = 1 px) and spatially down-sampled by a factor of 2. The final resolution of the measurement time series was thus 64 x 84 px (with each pixel being 19.4 x 19.4 µm) at 0.42 Hz. The concatenation of the pre-processed frames for all odours lead to the measurement matrix ∈ 5×7 with element 9,; being the observed value of the =0 pixel in the =0 frame.
Functional imaging (128 x 128 pixel, 5.02 µm/pixel) was performed for each odour 12.77 s with 94 Hz at a fixed z-position. To increase signal-to-noise ratio and reduce computational load we binned the raw data in time (120 frame bin, corresponding to a frame rate of 0.94 Hz). Then the odour induced activation was calculated as the increase in fluorescence − ' . ' was the mean fluorescence on the first 2 s after recording onset. Furthermore, the data was low pass filtered (Gaussian filter = 1.5 px ) and down-sampled by a factor of 2. The final resolution of the measurement matrix was thus 64 x 64 px at 0.94 Hz. To account for bleaching effects, the same procedure was applied to control measurements without any odour stimulation. Out of these we calculated for each pixel the mean bleaching time course and subtracted those from the pixels' odour responses time courses.
Before and after functional imaging we acquired anatomical images (512 x 512 pixel, 1.25 µm/px) at 3µm steps in z-direction. After functional imaging the mice were first killed by an overdose of urethane. Subsequently the OB was cleaned with artificial cerebro-spinal fluid

Glomerular response spectra
We segmented the functional image series into individual glomeruli and their activation course by means of regularized non-negative matrix factorization (rNMF) 36 . Such a factorization disaggregates the measurement matrix into components with a spatial signal distribution D and a common activation course D of the participating pixel. The factorization contained spatial smoothness (governed by parameter GH ) and sparseness regularization ( G; ) to promote a disaggregation into spatial distinct but possibly overlapping components, i.e. glomeruli arranged side by side with an overlapping signal distribution.
We decomposed IOS image series into 150 components with fixed smoothness regularization parameter GH = 2. The sparseness regularization parameter G; was adjusted such that spatial component correlation was just below 0.5 (see Soelter et al., 2014, for a justification of this value). In image analysis performed to extract MOR18-2 responses, we dismissed the usual non-negativity constraint on the components' activation courses to capture possible odour induced inhibitions. In contrast, in analysis to extract response spectra across the full dorsal OB we kept the non-negativity constraint as it increased the trial-to-trial correlation of response spectra. To remove non-glomerular components of the factorization, we only took components into account that exhibited a trial-to-trial correlation above 0.6 between response spectra of odour set repetitions within an animal 36 .
Components were manually assigned as MOR18-2 glomeruli if their pixel participation matched the GFP marked location of MOR18-2. Since the overall response strength varied strongly across animals we normed MOR18-2 responses to its Methyl propionate (MP) response in the same measurement ( OP*Q7 = / Q7 ). The final odour response was than computed as the median response of all measured animals, in which sufficient strong responses were observed ( Q7 > 0.2‰).
To evaluate similarity of glomerular response spectra we computed the pairwise correlation distance, O (eq. 1): with X + the response of glomerulus to odor , and 〈 X 〉 + the average response of glomerulus to all odors. Thus, highly correlated response spectra will have distance zero, whereas fully anticorrelated spectra will have distance 2.
Based on this distance we obtained a hierarchical clustering using the unweighted pair group method with arithmetic mean (UPGMA) approach and evaluated clusters at different thresholds. We manually set thresholds to obtain clusters with distinct gaps between them. Thereby we restricted clusters to contain at least glomeruli of three different animals. To account for varying overall response strength across animals we normalized all responses to unit length response spectra b+OH + = + /‖ ‖. by sampling from this spectrum from = 0 to = 4000 in intervals of . We fitted activation models m + = ( ) descriptor sets to match the observed odour activations + . We evaluated goodness of fit by calculating the coefficient of determination _ , eq. 2,

Physico-chemical characterization
over all odorants . Whereas this measure assessed how well the data could be described by a model, it does not reflect that sufficient complex models could perfectly fit any functional relationship of the data without the ability to generalize to unseen data. Therefore we also evaluated the coefficient of determination for unseen data in a bootstrap approach 38 . To this end, we fitted each model on 50 bootstrap samples 1 of the data and obtained predictions of odours m + not contained in the bootstrap sample. We then calculated the predictive power as in eq. 3, ( We evaluated the performance of univariate linear regressions, isotonic regression and unimodal isotonic regression 39 and multivariate Support Vector Regression (SVR) with Gaussian Kernel 40 . All models were obtained via the scikit-learn module using default parameters 41 .
Scikit-learn was also employed to obtain the Multidimensional Scaling (MDS) of odour space distances. Prior of fitting SVR and MDS models we normalized all descriptors to zero mean and unit variance with respect to our in-house library of about 1000 odours.
We also employed normalized descriptors to obtain glomerular position in chemical space. To this end we assigned each glomerulus its barycentre of activation, eq. 4, relative to the mean descriptor value ⟨ ⟩ + = L N ∑ + + for all odours in the measurement set. We then calculated the pairwise cosine distance of the barycentres (eq. 5), and obtained a hierarchical clustering thereof using the UPGMA approach.

Replicability
We wish to encourage further use of the dataset and methods developed here. Therefore, we provide all code and data that is necessary to replicate each figure in this contribution. Please refer to the "Replication guide" in the supplementary material for links and instructions.

Results
Our aim was to assess chemotopy in the dorsal olfactory bulb. We started by determining a detailed receptive range of the MOR18-2 glomerulus, from which we inferred a physicochemical model of the structure-activity-relationship for MOR18-2. The odorant profile and the physicochemical model enabled us to assess similarity of response profiles of glomeruli in the neighbourhood and evaluate the chemotopic order among glomeruli.

The molecular receptive range of MOR18-2
We measured the response of MOR18-2 glomerulus by Intrinsic Optical Signal (IOS) imaging of the dorsal olfactory bulb (dOB, Fig. 1a) under stimulation with monomolecular odorants in MOR18-2-IGITL mice 34 . The GFP label enabled direct localisation of the glomeruli on which MOR18-2 OSNs converged (Fig. 1b). Band-passing the signal reduced signal noise and separated the glomerular signal from the global unspecific signal 42 (Fig. 1c). Signals from individual glomeruli were then extracted using regularized Non-negative Matrix Factorization (rNMF)   In order to obtain an as complete as possible MRR of MOR18-2 glomeruli, we pursued an iterative screening approach, alternating between biological measurements and virtual screening 44,45 . In each iteration, we built computational models based the physicochemical properties of the biologically confirmed ligands. We used this model to screen odorant databases for potential ligands which we subsequently tested in vivo.
After discovering the first ligand (Ethyl acetate, compound 2 in Table ST1) by screening odours known to activate glomeruli in the dorsal bulb, we iteratively built hypotheses for further ligands based on the present set of revealed ligands. We acquired those hypotheses both empirically, by an intuitive notion of chemical similarity (e.g., chain length or bond saturation), and automated by virtue of virtual screening, using the computational models we built on the basis of physicochemical descriptors.
After several iterations of virtual/biological screening, we had obtained the response spectrum of the MOR18-2 glomerulus for 213 odorants from measurements in a total of 41 animals ( Figure 2a). The response spectrum of MOR18-2 showed a rather narrow range of strongly activating odours. The structural diversity of the MOL18-2's MRR is shown in Figure 2b, that depicts the most active ligands as well as some chemically similar but non-activating molecules. MOR18-2 was most sensitive to small aliphatic esters like Methyl Propionate, (compound 5), Ethyl Acetate (2), Propyl Acetate (4), Methyl Acetate (5), and the corresponding propionates (3,6,10). However, the most potent ligand we tested was Allyl Acetate (1), that has a double-bond in the terminal group. Interestingly, Isopropyl Acetate (9) yields relatively strong activity. It has a branched aliphatic group in the position of the double bond in Allyl Acetate.
The structural variations of Allyl Acetate (1) and Isopropyl acetate (9)    We also observed some small, but statistically significant negative responses (e.g. compounds (187) and (201)). This observation may indicate that the glomerulus innervated by MOR18-2 receptor neurons shows inhibitory responses to these compounds. Since we assume that the intrinsic signal is dominated by ORN input, this observation implies that MOR18-2 normally emits a certain base firing rate that is decreased by these compounds. However, we cannot rule out the possibility that bandpass filtering in pre-processing induced artefacts that appear as inhibitory responses, for example through a very strong responses of a closely neighbouring glomerulus.
Glomeruli are not always arranged in a single plane in the rodent olfactory bulb. In order to exclude the effect of glomerular overlay an imaging method with 3D-resolution was therefore necessary. As an additional validation that the extracted IOS signals reflect mainly the glomerular input activation to the MOR18-2 glomerulus, we performed high resolution 2-photon synapto-pHluorin (spH) imaging in mice expressing both spH in all OSNs as well as GFP in MOR18-2 OSNs. To this end we first performed a 3D-anatomical scan of the resting spH fluorescence in which we manually outlined glomerular contours (Fig. 3a). We then performed functional spH imaging (Fig. 3b) and extracted glomerular activation via rNMF (Fig. 3d).  The extracted spH activation of MOR18-2 glomeruli strongly affirmed the IOS derived ligand rank (Spearman's rank correlation 0.92, < 10 }~) . But in contrast to a rather continuous increasing response strength in the IOS, the spH derived spectrum had a pronounced separation of strong and weak activations (Fig. 3e). This might be attributed to either a steep non-linear gain of spH fluorescence, or to a saturated IOS responses for strong activations (or a combination of both). In any case, very high rank correlation of both spectra justify to rely on IOS measurements for subsequent analysis, which allowed us to record larger spatial extent and increased number of measured stimuli/repetitions than spH-imaging.

Physico-chemical receptive range
Having obtained the MRR of MOR18-2, we asked if it could be described by a quantitative ligand-based model, that is, if response strength could be explained by a physico-chemical model of the odorants. To this end we first obtained 1600 physico-chemical descriptors via eDragon 37 . These descriptors are subdivided into 19 blocks, ranging from simple scalar representations of molecules, like molecular weight or functional group counts, up to representations of the three-dimensional arrangement of properties (e.g. mass, charge) within the molecule. Additionally, we calculated the EVA descriptors in a range of spectral resolutions 46 . We used EVA descriptors previously to determine structure-activity relationships for olfactory receptors 23 . EVA descriptors have a free parameter, i.e. the width of the Gaussian kernel used for smoothing the vibrational line spectrum. We denote the particular variant used as EVAx , with x the kernel bandwidth in cm -1 (see methods for details).
We employed Support Vector Regression (SVR) models to explain and predict glomerular activation. Since we found that single descriptors are generally insufficient to describe structureactivity relationships in our data (see supplemental Results and Figure S3), we employed a multivariate approach. In a broad sense, SVR models estimate the response of the glomerulus to an odorant by forming a weighted sum of the responses to similar odorants, where similarity is determined with the underlying descriptor space. Therefore, for SVR models to become predictive, odorants evoking similar glomerular response must be proximal along some axis of the descriptor space. A multidimensional scaling (MDS) projection of odorants in the combined space of all eDragon and EVA5 descriptors showed that activating odorants clearly cluster together (Figure 4a). Consequently, the SVR model obtained in the combined eDragon/EVA5 space was able to capture the glomerular response to odorants in the training set ( Fig. 4b), as well as generalise to odorants in the test set, which were not involved in model construction in any way (Fig. 4c). We evaluated which physicochemical descriptor sets would yield the highest predictive power in SVR models (Fig. 4d). We compared the full set of eDragon descriptors (eDRAGON), the logical subsets within them (see supplementary table ST2 for abbreviations), two subsets previously proposed to describe odour distances (HADDAD 21 , SAITO 20 and the EVA descriptor at resolutions of 1, 5, 10, and 50 cm -1 (EVAx). We found that the full set of eDragon descriptors yielded the best performance (followed by the 2DAUTO (2D-autocorrelations) and the odour optimized HADDAD subsets. Notably, while the EVA descriptors did not outperform the HAD-DAD and SAITO descriptor subsets, they still yielded better predictions than most of the other eDragon subsets. We then tested whether predictive performance could further be enhanced by combining the eDragon descriptor set with the EVA descriptors (eDR-EVAx, Fig. 4e). And indeed, the combination of both sets improved performance by up to about 30 %, yielding models of MOR18-2 activation with a predictive power _ ≈ 0.4. Thus we conclude that, despite its high dimensionality (2400 descriptors), proximity in the combined eDragon-EVA descriptor space provided the best estimate of chemical similarity of odorants with respect to their activation of MOR18-2.

Tunotopic embedding
Since our imaging approach covered a considerable field of view across the dOB, it also included simultaneous recordings of glomeruli in the neighbourhood of MOR18-2, allowing to compare its MRR to other MRRs in the local glomerular ensemble. In general, the main class of MOR18-2 ligands (i.e., small esters) activated several glomeruli across the dOB (Fig. 1b and Supplementary Figure S1). Activity elicited by MOR18-2 ligands was therefore not restricted to a clearly defined sub-region. Likewise, the spatial location of the MOR18-2 glomeruli was not designated by any odorants which were solely activating this area. This raised the question about how the molecular response profile of the MOR18-2 glomerulus related to the tuning of its neighbours -that is, its tunotopic embedding in the dOB.
To investigate this tunotopic embedding we compiled a set of 48 odorants that we used to systematically probe the MRR of the glomeruli surrounding the MOR18-2 glomerulus (Table   ST1, last column). This set contained most of the ligands that elicited a significant positive response in MOR18-2, along with some of their structural variations, therefore representing the previously determined MRR of MOR18-2. We measured the dOB IOS response in four GFP and three SpH mice and extracted glomerular response spectra (see methods). In the SpH mice, we validated the glomerular segmentation of the IOS via anatomical glomeruli outlines obtained in a 3D-scan of the resting SpH fluorescence (see 36 ). The GFP marker allowed us to identify the MOR18-2 glomeruli.
To investigate tunotopic relations between the extracted glomeruli, we calculated pairwise response similarity between all extracted glomeruli in all seven animals. Response similarity was assessed as the correlation distance O of the odour response spectra. Based on this distance we then obtained a hierarchical clustering of glomeruli (Fig. 5a). The spectra of the GFPlabelled MOR18-2 glomeruli displayed high mutual similarity, and they were all contained in a cluster with a low similarity threshold (red cluster in Fig. 5), suggesting that this approach is suitable to identify glomeruli corresponding to the same olfactory receptor across individuals.
We observed several other low-threshold clusters that indicate a putative match of glomeruli across animals, indicated by the yellow, green, turquoise and blue clusters in Fig. 5. The glomeruli corresponding to those clusters were located at stereotypic positions across animals ( Fig. 5b and Supplementary Figure S2). Their individual response spectra were highly overlapping (Fig. 5c). Histograms of pairwise similarity values of these glomeruli to all other glomeruli clearly showed that the putatively matching glomeruli shared high similarity to each other, while the similarity to non-members of the clusters was consistently lower, at least for the red, yellow and green clusters (Fig. 5d). For the turquoise and the blue clusters, the observed response spectra were not as clearly separated from the other glomeruli as for the red, green and yellow clusters. Nevertheless, a clear separation was possible when assessing the similarity to the cluster centroid, instead of a full pairwise comparison (Fig. 5e). In the case of the green cluster in Fig. 5, the cluster spanned only a subset of the animals -that is, the corresponding glomeruli could not be extracted in all 7 animals -probably, because it was located outside the recorded area in some of the animals. Nonetheless, we conclude that glomeruli with high mutual similarity likely correspond to glomeruli driven by the same class of olfactory receptor neurons, and that they can be identified by the similarity of their functional response spectra through hierarchical similarity clustering. We next investigated the similarity response profiles of those putatively identified glomeruli to the MOR18-2 glomerulus. To this end we determined the correlation distance † ‡ˆL‰}_ = ( † ‡ˆL‰}_ , 1 ) between the spectrum of a putative glomerular cluster 1 to the spectrum QŠ‹L‰}_ of the cluster corresponding to MOR18-2.
In Fig. 6a, clusters are highlighted if they exhibited a response correlation QŠ‹L‰}_ > 0.2. The hierarchical clustering procedure grouped four of those clusters together with the MOR18-2 cluster (coloured blue to violet in Fig. 6a), forming a 'meta-cluster'. Closer inspection of their response spectra revealed that these glomeruli share an activation by propionates (Fig. 6c).
These glomeruli did not only share similar activation but were also spatially proximal to each other on the dOB surface, forming a patchy domain around the MOR18-2 glomeruli (Fig. 6b).
Besides this local domain of glomeruli with mutually similar response spectra, MOR18-2's spectrum was also correlated to other putative glomeruli located lateral-posterior to MOR18-2 glomeruli. However, each of those cluster was most similar to glomeruli uncorrelated to MOR18-2. That is, although glomeruli with somewhat similar response profiles to MOR18-2 were located at distant locations, only those in proximity to MOR18-2 formed a meta-cluster of high mutual similarity. This clustered spatial arrangement was found in all animals under consideration.
We then asked whether the relative arrangement of such meta-clusters of spatially grouped glomeruli also followed a stereotypic rule that was conserved across animals. To this end, we chose a high similarity threshold to form clusters of glomeruli, obtaining essentially two 'superclusters' of glomeruli (Fig. 6d). The members of one supercluster were located predominantly in the lateral-posterior part of the dOB (blue glomeruli in Fig. 6e), while the other supercluster's members were located mostly in the medial-anterior part (coloured red in Fig. 6d and e). The members of the MOR18-2 group (coloured green in Fig. 6e) are located in a region where those two 'superclusters' overlap. This pronounced bipartite arrangement of glomerulus superclusters in the dOB was observed in all 7 animals that we analysed. In summary, we found that glomeruli were arranged at stereotypical positions regarding their response similarity to MOR18-2. Moreover, their spatial arrangement in the dOB correlated with the chemical response profile of the corresponding olfactory receptors. This holds on a local level, as a group of glomeruli with high overlap in their response spectra were arranged in a confined, but patchy, spatial domain around MOR18-2. Furthermore, we found that, on a larger scale, glomeruli in the dOB can be grouped in two 'superclusters' according to their response spectra, separated in a lateral-posterior, and a medial-anterior domain of response similarity.

Chemotopic embedding
The tunotopic embedding of MOR18-2 that we observed suggested that there would exist, along some dimension of chemical space, also a chemotopic embedding, that would group glomeruli with a similar physico-chemical response profile profile together. Testing for chemotopy is, however, notoriously difficult since there is yet no clear-cut low-dimensional description of odorant space that would allow to define a sensory topology of the like that is found in visual or somatosensory systems. Therefore, investigating a potential chemotopic embedding required defining a chemical space (i.e., chemical similarity) in which odorants are arranged, and to define the position of glomeruli within this space.
As chemical space we choose the combined eDR-EVA5 space since it proved best to characterize the MRR of MOR18-2 (see section MOR18-2 physico-chemical receptive range). Furthermore, we placed the origin of the coordinate system at the mean descriptor value ⟨ ⟩ + of our odour set. The location of a glomerulus within this space was defined by its barycentre (i.e., the centre of mass) D with regard to its activation spectrum (see eq. 4 in the Methods section). Hence, a glomerulus responding with equal activity to all odorants would be positioned at the origin of the chemical space.
We then obtained the relative positioning of glomeruli to each other by calculating their cosine distance •OE• = D ⋅ * , which is a measure of the angle between the glomerular positions (see eq. 5 in the Methods section). In other words, glomeruli were considered similar, if their barycentres of activation were located in a similar direction from the origin.
Based on these pairwise distances we again obtained a hierarchical clustering (Fig. 7a). All MOR18-2 glomeruli were grouped together in a low threshold cluster (violet cluster in Fig. 7a & b), similar as with tunotopic clustering (cf. red cluster in Fig. 6a & b). We also found, again, that glomeruli most similar to MOR18-2 were predominantly located in its immediate vicinity, or slightly lateral-posterior to its location. However, in contrast to the tunotopic case, the MOR18-2 glomeruli were not embedded in a small scale local chemotopic domain, but only in a meta-cluster of lateral-posterior located glomeruli (violet-blueish cluster in Fig. 7c & d). This meta-cluster was contrasted by an adjacent meta-cluster of medial-anterior glomeruli (orange and yellow cluster in Fig. 7c & d). But on the topmost level of clustering more distant glomeruli were widely dispersed in both domains (grey glomeruli in Fig. 7c & d). In summary, our approach could reveal a clear chemotopic arrangement of glomeruli with high similarity to MOR18-2. However, the observed chemotopic structure was less obvious than when assessing response similarity on a tunotopic basis.

Discussion
In our study, we profiled in detail the MRR of MOR18-2 glomeruli and their embedding in the glomerular code. MOR18-2 glomeruli mainly responded to small esters and locally embedded in a patchy domain of glomeruli sensitive to short chain propionates. Furthermore, we could determine a physico-chemical receptive range of MOR18-2 glomeruli within the multidimensional eDR-EVA5 space via SVR models. We also observed a chemotopic representation of this space within the dorsal OB, albeit less obvious than the tunotopic arrangement.

MOR18-2 receptive range
We aimed to obtain the MRR of MOR18-2 as completely as possible, but it is of course infeasible to screen all existent odours. Our informed search for candidate ligands was based onpredictive physico-chemical models, but also (to a lesser amount) randomly selected search candidates. In spite of this directed search, it is unlikely that our set of ligands exhaustively represented the full response repertoire of MOR18-2. Indeed, our predictive SVR model, as ascertained on bootstrap samples, suggests further possible ligands, which could be evaluated in future screening campaigns. To our knowledge this study makes MOR18-2 one of the best explored olfactory receptors 47 . Based on the data published along with this report, it will be possible to construct new predictive models, potentially including any ligands discovered later.
Recently, acetic and propionic acid have been reported as ligands of MOR18-2 receptors expressed in the kidney 48 . In contrast, our results rather identified the corresponding esters as ligands. Although their study did not contain any of the ligands revealed in our study, there could be systemic reasons for such a deviation. For example, we presented odour stimuli in the gaseous phase, taken from the headspace of a solution, while Pluznick et al. applied the liquid phase. Since acids tend to have much lower vapour pressure, the effective concentration that reached the receptors might have been considerably lower for acids than for the corresponding esters. Another potential reason for this discrepancy is the potential conversion of esters into the corresponding alcohols and acids by enzymes in the mucus 49 .
Assuming that esters are generally converted to their corresponding acids shines an interesting light on our results. The best 12 ligands we identified are all esters of either acetic or propionic acid. Hence, all of their responses could be explained by the receptor responding to either of those acids. Asserting whether this is indeed the case would require administering esterase blockers in the mucosa prior to the experiment, which is beyond the scope of the current study.
As our aim was to obtain an as complete MRR of MOR18-2 as possible, we did not obtain doseresponse curves for odorants, since obtaining those curves would have dramatically reduced the amount of odorants to be measured. Instead, all odorants have been used at very high concentrations -liquid odorants in highest purity available and solid ones in saturated solutions. We are therefore unable to make any statement of concentration effects on the MRR.

Physico-Chemical Receptive Range
In this study we did not only determine the MRR of MOR18-2 but also translated it via SVR models to a physico-chemical receptive range in the eDR-EVA5 space. Despite bootstrap validation of the model and the success of a very similar approach for Drosophila ORs 23 , an ultimate confirmation by successful biological screening of further model predictions is still missing. Nonetheless our careful model validation ascertains that our obtained physico-chemical range is at least a rough estimate.
We obtained the best results in describing the physico-chemical receptive range by combining features based on molecular vibrations (EVA) with features containing descriptions of the molecular shape (eDragon). Note that this provides neither evidence for nor against the highly disputed theory of vibrational olfaction 50 . In our model, molecular vibrations are considered as a kind of "fingerprint" of the molecular graph, rather than a mechanistic interpretation of molecule-receptor interaction.

Tunotopic domains
Our results provided evidence for tunotopic domains at different scales: two global domains,  51 . It is also in line with the MOR18-2 glomeruli being located close to the border of class 1 and class 2 domains 52 . In contrast, the local tunotopic clusters we identified diverge from those identified in other studies, probably because of the different odorant sets that were used to identify tunotopic groups of glomeruli.
Interestingly, despite being a class 1 receptor, MOR18-2 is tunotopically embedded in the lateral-posterior domain which is presumably composed to a large extent of class 2 receptors.
This could indicate that, albeit its physiological separation 52 the functional transition between receptor class domains is rather continuous.
Mori and colleagues 14,51 also observed a compartmentalization of the dOB into local response clusters. They attributed long chain ester responses to cluster A (within class 1 receptor domain) and cluster B (within class 2 receptor domain) but did not attribute any cluster to small esters. However, their cluster assignments leave undesignated space especially at the border of domain 1 and 2 where we observed the propionate responding cluster of MOR18-2. This area is potentially populated with TAAR glomeruli 53 .
Taken together, our findings confirm the hypothesis of a tunotopic layout of the olfactory bulb.
This was disputed by 17 who observed "that nearby glomeruli were almost as diverse in their odorant sensitivity as distant ones". In their study, they evaluated the dependence of glomerular spectrum similarity with regard to their spatial distance and did not find strong effects. A potential reason for the different result could be that a highly diverse odour set such as the one used on that study is not suited to resolve the similarity in MRRs among local response clusters. For example, the local response cluster of MOR18-2 would not emerge using an odour set containing only acetates, but lacking the corresponding propionates (which differ only by one additional Carbon atom).

Chemotopic domains
Besides tunotopic domains, we did also observe some locally confined cluster of glomeruli with chemical similar response spectra. This chemotopic organization was less pronounced than the tunotopic organization, for several potential reasons. First, there is no canonical definition of chemical similarity. Our definition by the eDR-EVA5 descriptor space is suitable to describe the MRR of MOR18-2, but it might still be too imprecise to reveal physico-chemical structure-activity relationships in local chemotopic domains. Second, since our odour set was constructed to contain MOR18-2 ligands, it naturally would not cover the full ligand spectrum of other receptors. While we made this choice deliberately in order to resolve the local domain (in a chemical sense) of MOR18-2, our set of odorants will not be optimal to describe the MRRs of its neighbours. Since we represent the chemical receptive field of all glomeruli as a weighted combination of their response to the measured odorants, a full comparison of those receptive fields would require the key odorants of neighbouring glomeruli also being present in the set.
This does not so much pose a problem when assessing tunotopic distance, since the response spectrum correlations are only influenced by the identity, but not the chemical properties of the ligands. In contrast, chemotopic distance is all about the chemical properties of the shared and unshared ligands. A complete assessment of chemical similarity between MRRs can, in theory, only be assessed if both MRRs are represented equally well a number of ligands that sufficiently represent the region of chemical space in which those MRRs are located. Future research in this area might focus on identifying the MRR of neighbouring glomeruli to a much higher extent than in the present study, thus enabling better resolution of chemical similarity between MRRs.
A rigorous quantitave assessment of chemotopy is further complicated by the ambiguity of glomerulus locations. Most animals tested had multiple instances of the MOR18-2 glomerulus that were in the same area, but not immediately adjacent. In that case, it is impossible to provide a distance metric that appropriately captures chemical distances between MOR18-2 glomeruli and any other identified glomerulus. We considered using the centre of gravity of MOR18-2 locations as a reference, or running separate assessments for all instances of MOR18-2 glomeruli. Both approaches will result in non-zero distance between glomerulus instances that should be treated as identical for the purpose of mapping. As this violates the definition of a distance metric, any quantitative analysis based on such an approximation would be handwaivy.
To conclude, we investigated the olfactory bulb's topographic layout with regard to an extensively characterized MRR of a single glomerulus. This provided further evidence for a functional topographic layout according to glomerular response profile similarity. Nonetheless also this approach still suffered from the incomplete representation of all but one MRR. Therefore, it once more emphasises the demand of a full characterization of olfactory MRRs.

Additional information Funding and Acknowledgements
This works has been enabled by the Priority Programme SPP1392 "Integrative Analysis of Ol