Transfer Function-Based Characterization of the Honey Bee Olfactory System: From Biology to Electronic Circuits

We exemplify an interdisciplinary approach wherein a mesoscopic-scale functional model of a biological system is derived from time-series recordings, yielding transfer functions that can be used to design analog electronic circuits. Namely, sensory processing in the honey bee, a universal model for studying olfaction, is considered. Existing studies have focused on its antennal lobe, wherein only the responses of its functional units, known as glomeruli, have been accessible. Here, high temporal resolution calcium imaging is deployed to track the dynamics of odor-evoked activity beyond this processing stage. The responses in the somata outside of the antennal lobe are recorded, showing for the first time how the glomerular signals are transformed before entering the higher brain centers. A transfer function approach is applied to capture as a “gray box model” the remarkably heterogeneous signal transformations between odor input and glomerular response, and between glomerular signals and somata activity. The somata are tentatively mapped to the glomeruli via Granger causality, while machine learning classification and clustering allow grouping common properties regarding response amplitudes and temporal profiles. The obtained low-order transfer functions display time- and frequency-domain input-output properties closely similar to the biological system. Because transfer functions have universal applicability, once they have been determined, it is readily possible to design corresponding analog electronic circuits, with possible future applications in sensor signal conditioning. To exemplify this, examples based on resistor-capacitor (RC) networks and operational amplifiers are physically built and confirmed to generate responses highly correlated to the initial biological recordings.


I. INTRODUCTION
This proof-of-concept work aims to introduce and exemplify an analysis and design approach based on extracting mesoscopic-scale calcium imaging time series from neurophysiological experiments and deriving transfer functions The associate editor coordinating the review of this manuscript and approving it for publication was Wei Wei. from them. In turn, these allow designing physicallyrealizable analog electronic circuits mimicking the responses of neural ensembles, herein focusing on the paradigmatic case of olfaction.
Olfaction is a neurobiological process involving a complex sequence of signal transformations, starting from plumes of molecules in the environment having a spatiotemporally varying concentration. These are chemically sensed and undergo multiple representations, eventually leading to odor recognition and informing an organism's decision-making. Owing to its fundamental nature and evolutionary importance, olfaction has been the subject of intense experimental and theoretical analyses, which have led to manifold models across species. These are diversified by their level of physiological detail at the micro-, meso-and macroscopic scales, and primarily feature a feed-forward topology [1]- [8]. So far, the olfactory system has been extensively investigated primarily in the fruit fly (Drosophila), the microcircuits of which have been mapped down to considerable detail [9], [10]. For the case of the honey bee, herein taken as a ''toy model'' example towards a possibly general approach, it is established that the stimulus parameters are firstly encoded within specialized neural structures known as the antennal lobes (ALs) according to stereotyped activation patterns. These are subsequently forwarded to higher-order brain centers, such as the so-called mushroom bodies (MBs) and lateral horns (LH), where they are decoded and evaluated, potentially triggering a behavioral response ( Fig. 1) [11].
At first, in the olfactory receptors on the antenna surface, molecular concentrations are transduced via a ligand-receptor binding process into the electrical activity of olfactory receptor neurons (ORNs) [12]. This step can be represented as an initial low-pass filter, which is followed by a differentiating linear filter that transforms the membrane potentials generated by receptor activation into trains of spikes, in turn conveyed to the AL [13]. The AL is formed by neuronal ensembles known as glomeruli, each of which receives its input exclusively from a single type of ORN. Within the glomeruli, the afferent signals from the ORN axon terminals elicit postsynaptic potentials in the dendritic fibers of projection neurons (PNs). In this way, the original spike sequence is converted into fluctuations of local membrane potentials, which are modulated by the local neurons (LNs), that can be excitatory (eLNs) or inhibitory (iLNs) [14]. The excitatory modulation by eLNs is fast and accounts for the spiky initial phase of the glomerular PN response.
Contrariwise, the inhibitory action of the iLNs occurs slowly and is responsible for attenuating the response observed over time [15], [16].
Throughout these steps, the signals are reshaped in both their amplitudes and temporal features [14], [17]- [19], which together encode odor identity and intensity [20]- [22]. Finally, the signals from the PN dendrites are transformed back into spike trains [23] and propagate towards higherorder brain centers. Directly measuring a traveling signal simultaneously at different locations across a multistage neuronal circuit poses considerable instrumental difficulties, and, consequently, relatively little is known on how exactly the input signal is transformed along this path [11]. Here, we attempt to address this question at two levels. First, leveraging the increased temporal resolution of two-photon microscopy, we simultaneously probe the glomeruli within the AL and the somata outside it. Second, we model the resulting multivariate data with an approach centered around standard systems engineering notions, namely, using transfer functions to capture the input-output relationship across each processing stage.
By construction, suitable transfer functions can encapsulate an arbitrary device response. Therefore, while somewhat removed from biology, they are inherently more generally applicable than ad-hoc models such as those based on spiking networks. Transfer functions represent inputoutput relationships in the complex Laplace domain through the variable s [24], without needing to account in detail for the system internal variables. Consequently, based on potentially heterogeneous sets of experimental input-output observations, they allow extrapolating simulated responses to idealized impulse and step inputs, which can be compared directly. In brief, when restricted to the imaginary axis, i.e., s = 2π √ −1f , the abstract variable s allows to compactly represent the frequency response to purely sinusoidal tones (i.e., at frequency f ), thus accounting for the complex Fourier domain. More in general, the roots of the numerator and denominator polynomials (known, respectively, as the system zeros and poles) determine with considerable generality the characteristics of the time-and frequency-domain responses. For example, complex conjugate poles with low damping give rise to oscillatory dynamics in the time-domain response, and to resonance peaks in the frequency response.
In the present study, transfer functions are particularly suited for the purpose because they provide an external, though partially interpretable, model of the system in terms of how it responds to inputs. Transfer functions combine a coarse-grained theoretical understanding of the underlying process, with experimental data to determine the free parameters. As such, they make be considered a kind of ''gray box'' model, representing an intermediate scenario between a non-interpretable purely data-driven approach, such as a neural network learning a relationship without making prior assumptions about its form (black box), and a complete, detailed mechanistic account derived from fundamental principles and physical laws (white box). This approach is FIGURE 2. Experimental setup. a) An olfactometer produced odor pulses while the neural responses were synchronously recorded via a two-photon microscope from a bee mounted on the imaging stage. b) Independent photomultipliers (PMT) in two detection channels set up so as to allow simultaneously recording the calcium signals from the glomeruli and somata despite their substantial intensity difference.
Transfer functions do not require any knowledge of the system internal variables, which, indeed, in the case of a mesoscopic neural structure such as the antennal lobe are not fully accessible. Despite their versatility, transfer functions have thus far been deployed to neuroscientific problems only sparingly, plausibly because of their inability to represent learning or adaptiveness, only providing a filter-based simplification of a system. Nevertheless, their empirical usefulness has been confirmed in capturing the responses of single neurons [27], functional interdependencies between brain regions, and neurovascular couplings in functional MRI [28]- [30]. Here, we apply transfer functions to decipher the olfactory transcoding process into a representation suitable for informing the engineering of a signal processing model or even another physical device [31], [32]. In doing so, we present an example of a systematic approach that, starting from biological recordings, eventually yields analog electronic circuits that are compact and straightforward to realize physically. We therefore posit that, in addition to their usefulness as models of neural activity, transfer functions are helpful towards designing bio-inspired sensing devices mimicking a subsystem of a natural organism, albeit with some important limitations.

A. SPECIMEN PREPARATION FOR in vivo CALCIUM IMAGING
Fourteen honey bees (Apis mellifera) were exposed to CO 2 for 30 s, then fixed onto a custom-made imaging stage (Fig. 2a) using soft dental wax. After that, a small rectangular window was cut into the cuticula, and the glands and trachea covering the AL were gently moved aside. Subsequently, fura-2-dextran, a calcium-sensitive fluorescent dye (type F3029, Thermo-Fisher Scientific, Inc., Waltham MA), dissolved in distilled water, was injected into the antenna-cerebral tracts, posterolaterally to the α-lobe, using a pulled glass capillary. After injection, the excised cuticula was fixed back in its original position using n-eicosane (Sigma-Aldrich, Saint Louis MO). The bees were stored in a dark, cool, and humid place for about 20 hours to allow dye diffusion into the AL.
The invasiveness of the preparation and experimental procedures was kept to a minimum. Furthermore, the integrity of the nervous system was repeatedly tested behaviorally via the proboscis extension reflex after contact stimulation of the antenna with a sugar solution, largely ruling out confounding effects of damage [33]. During recording, the stability of the neuronal responses as visible to calcium imaging was monitored by an experienced operator, and no instances of response loss or abnormal drifting were detected.
All procedures were completed according to approved institutional protocols of the University of Trento (Trento, Italy). In accordance with applicable laws, no project-specific institutional or national approval was needed for experimental work on insects.
Before the imaging session, the antennas were blocked with a drop of n-eicosane on the pedicel, leaving the flagellum free to move. The cuticular window, the trachea, and the glands were moved away from the antennal lobe region. A silicone adhesive was used to cover the brain, and a rectangular plastic foil was attached frontally to the window in order to separate the antennas from the immersion water required by the objective lens.

B. TWO-PHOTON MICROSCOPY
The two-photon microscope (Ultima IV, Bruker, Inc., Billerica MA) was illuminated by a Ti:Sa laser (Mai Tai Deep See HP, Spectra-Physics, Inc., Milpitas CA), which was tuned to 780 nm for fura-2 excitation. All images were acquired using a water immersion objective (10×, NA 0.3; Olympus, Inc., Tokyo, Japan). The absolute fluorescence signal from the glomeruli of projection neurons (hereafter, for brevity, somata) is about 10-fold more intense than the one from the somata , due to different partial volume effects between the larger and denser glomeruli, and the smaller and sparser somata. Therefore, two separate detection channels were prepared to record from both with sufficient signalto-noise ratio (SNR). Namely, the received light was split into two branches, illuminating different photomultipliers, VOLUME 10, 2022 by means of a dichroic mirror. The lower band was filtered at 525 ± 35 nm, close to the fluorescence peak of fura-2, and used for imaging the glomeruli. The upper band >560 nm, corresponding to the tail of the fluorescence spectrum, was used without additional filtering to image the somata. This approach maximized the dynamic range without saturating the photomultipliers, and the resulting setup is visible in Fig. 2b. The laser power was set to ≈ 10 mW after the objective as a compromise between maximizing the SNR and limiting photodamage.
The image acquisition was synchronized to the stimulation protocol and performed at a frame rate of 9.7 fps. The resulting image of 128×128 pixels with a digital zoom factor of 3.8 covered a field of view of 280 µm. The fluorescence intensity was recorded with a depth of 13 bits. In addition to the functional images, for supporting the morphological identification of the glomeruli a z-stack of the antennal lobe was acquired at a spatial resolution of 512 × 512 pixels using a layer interval of 2 µm.

C. OLFACTORY STIMULATION
Stimulation was performed through a custom-built olfactometer [34], allowing the delivery of 8 different odorants in an air flow having a velocity of ≈ 4 m/s, corresponding to a gentle breeze. A pseudo-random sequence with millisecondprecision timing was used, and each stimulus lasted 3 s followed by a 12 s recovery period. The sequence was repeated over 10 trials. Multiple odorants were necessary to maximize ecological validity and reduce adaptation. All of them were diluted in mineral oil. For 3-hexanol, 1-heptanol, 1-hexanol, citral, 2-heptanone, and isoamyl acetate, the dilution was 1:200 (vol/vol). For benzaldehyde, the dilution was 1:50, and for geosmin, it was 10 −6 . These diversified concentrations were chosen specifically to study the bee olfactory system under naturalistic conditions, that is ensuring receptor activation while avoiding saturation, leading to the highest dynamic range and odor-specificity in the glomerular code.

D. SIGNAL POST-PROCESSING
Data post-processing and analysis were performed via custom scripts written in the MATLAB language (Mathworks, Inc., Natick MA). For each bee, the fluorescence signals were spatially averaged over and recorded from the glomeruli, identified using the AL atlas [35], and the neuronal somata, delineated using regions-of-interest drawn by an experienced operator. From these raw traces, the relative change of fluorescence during stimulation was calculated, and expressed as F/F = −[F(t) − F b ]/F b × 100, where F b denotes the average fluorescence signal intensity over the 3 s pre-stimulus period. Finally, for each glomerulus or soma, F/F was averaged over all trials. To identify the glomeruli having the highest variance during stimulation, a principal component analysis (PCA) was performed, taking the pixels as variables and the frames as observations [36].

E. RESPONSE CLASSIFICATION
The individual responses from the glomeruli and somata were separately classified as activation, inhibition, or absent by means of a machine learning approach. For this purpose, a range of classifiers and settings were empirically considered, and the highest-performing configurations according to 10-fold cross-validation were chosen. For the glomerular responses, a one-vs-one support vector machine (SVM) classifier with quadratic kernel and unitary box constraint was selected. For the somata responses, a k-nearest neighbors algorithm (kNN) with k = 10 neighbors, correlation distance as a metric and equal weighting, performed best. In both cases, data standardization was applied. In principle, SVM could perform better than kNN at this task, but it was observed that its generalization performance was lower for the somata, as indicated by cross-validation. This difference could be ascribed to the lower signal-to-noise ratio for the somata traces, a finding which is well in line with other studies that have identified superior kNN over SVM performance in noisy acoustic datasets [37], [38]. The underlying training sets are provided as Supplementary Materials.
The models were trained based on response classifications performed by an experienced operator. To adequately cover the stimulus and the immediate post-stimulus phase, over which some elicited activity persisted, 88 time-points (≈ 9.1 s), starting from 2 frames before stimulus onset, were entered in the analysis as predictors. Overall, 616 manuallyclassified responses were entered for training the model for the glomeruli, whereas for the somata, 1074 responses were available. As visible in Fig. 3, the resulting accuracy, intended as agreement with human classification, was 99% for the glomeruli and 83% for the somata responses.

F. TRANSFER FUNCTION ANALYSIS
To determine the transfer functions, and hence, to reveal the underlying heterogeneity of temporal dynamics, a further step of clustering was applied, this time within each response class (activation or inhibition) and location (glomeruli or somata). Following extensive consensus work (omitted for brevity), based on joint consideration of the Silhouette method [39], alongside aspects of the Calinski-Harabasz [40], and Davies-Bouldin [41] methods, a fixed number of k = 4 clusters was predetermined for all combinations.
This setting was arrived at as follows. Firstly, the Silhouette value S(k) was calculated while assuming k = 2 . . . Based on these values, k = 2 would be chosen as sufficient. However, as a second step, the corresponding average time-courses for the most tonic and most phasic responses were also visually compared ( Supplementary  Fig. 1). This revealed beyond doubt that the separation between the two increased substantially between k = 2 and k = 4, and considerably less visibly for higher numbers of clusters. Based on this observation, k = 4 was eventually chosen as a compromise between emphasizing the variability of the responses and limiting the number of clusters. None of the results presented hereafter are critically dependent on this choice; the transformation between tonic and phasic and viceversa would also be visible with k = 2, albeit less markedly (data not shown). This clustering step, therefore, should be intended as a convenient means of partitioning a continuum into discrete categories, rather than as aiming to uncover responses that are in themselves separated into groups.
The trial-averaged time series, each consisting of a central subset of 56 peri-simulus points (5.8 s span, chosen empirically), were therefore entered into k-means clustering based on squared Euclidean distance, with k = 4. The cluster sizes are reported in the Results section.
To determine the transfer function parameters within each cluster, the averaged response was considered. Furthermore, to reduce the dimensionality of the analysis, the responses from all odors administrated were pooled together, because no consistent odor-related differences between the responses could be decoded in the preliminary analyses. Although each odorant was administered at a single concentration, pooling the data from the different ORN featuring strongly differing sensitivities to the individual odors leads to response signals with very different amplitudes. Transfer function determination was performed separately for the two processing stages. Firstly, from the input odor time course to the glomerular response. Second, from the time courses of the glomeruli to those of causally-related somata, as detailed below.
To determine the input time series for the glomerular responses, the temporal profiles of the odorant concentrations were measured through photoionization detection (type miniPID 200B, Aurora Scientific, Inc., Ontario, Canada), except for 1-heptanol, citral, geosmin, and 2-heptanone, which were undetectable. The profiles were largely overlapping, indicating that the diffusion dynamics are primarily determined by physical parameters of the setup, such as gas flow rate. Therefore, using a stereotyped profile was acceptable, and preferable owing to the absence of noise. Consequently, all analyses assumed a stereotyped pulse, having exponential rise and decay respectively according to 1 − exp(−10t) and exp(−2t) (Fig. 4).
Transfer function estimation was performed assuming no delay, 3 poles, and up to 2 zeros. This parsimonious configuration yielded a good fit across all responses, as demonstrated below by the high correlation between the realized electronic circuits and the biological recordings. The initial conditions were automatically searched for using an empirical approach minimizing prediction error. A least-squares non-linear optimization for prediction error minimization was performed while enforcing stability and setting λ = 10 −8 during parameter regularization. The underlying input-output data and complete parameters obtained by the estimation process are provided as Supplementary Materials.
The fit quality was represented based on the Normalized Root Mean Squared Error (NRMSE) as Q = 1 − NRMSE. To summarize the transformations performed based on the step response, the static and maximum gains S 0 , S max were extracted, from which the relative overshoot was determined as σ = (S max − S 0 )/S 0 . Furthermore, the rise time τ rise was determined as the lag at which the normalized response reaches 1 − 1/e. It should be underlined that, by its nature, calcium imaging only offers an indirect measure of neural activity, and the amplitudes of the recorded responses are affected by numerous factors, including the possibly inhomogeneous level of dye uptake and, as elaborated in the Discussion section, the relationship of the signal to pre-and post-synaptic potentials across the glomeruli and somata. Therefore, the transfer functions empirically represent the transformation of the signal time-courses, but the corresponding static gains do not have a straightforward biological interpretation. Accordingly, in the below, all responses are graphed in unitless form, in contrast with electrophysiological recordings that offer calibrated current and voltage readings, and are necessary to properly determine the neural transfer functions [27], [42].

G. GRANGER CAUSALITY ANALYSIS
To date, it is technically not feasible to image in-vivo the structural connections between the glomeruli and somata. VOLUME 10, 2022 Therefore, we resorted to determining the functional interdependencies between them, informing at a statistical level which pairs to consider for determining the transfer functions between these two processing stations. Synchronization measures such as cross-correlation can be influenced by numerous confounding factors as exemplified in chaotic oscillators [43]. Thus, to obtain a more robust estimate also with respect to non-neural signal sources and other instrumental fluctuations, we applied Granger causality. This is a statistical measure deemed highly reliable in analyzing diverse neuroscientific datasets including EEG and functional MRI signals [44], [45].
Here, this analysis was systematically performed on the time series from each glomerulus-soma pair, namely, assuming each one as a possible driver and target within the individual odor response maps for all individual bees and odor stimuli. Each time series entered in the causality analysis consisted of the concatenation of 10 trials, considered separately for each odor, spanning 4 s (3 s during stimulation, 1 s post-stimulus). This approach was chosen over entering the entire time series in order to emphasize the interdependencies between the responses over the unrelated, stochastic fluctuations. To further reduce contamination and remove instrumental fluctuations together with other baseline instabilities, the signals were detrended and high-pass filtered using a second-order Butterworth filter at 0.1 Hz.
For each pair comprising a driver D and a target T, a linear regression of the present state of T, y T (t), was performed twice: once only given its past state Y − T (t) and yielding the prediction error e T|T (t) = y T (t) − E[y T (t)|Y − T (t)], and again, given the past states of both D and T, [46]. The corresponding prediction error variances, λ T|T = E[e 2 T|T ] and λ T|T,D = E[e 2 T|T,D ], together give the Granger causality, with Here, the past states were empirically approximated as To test for statistical significance, 100 sets of multivariate surrogate time series were generated using the Iterative Amplitude Adjusted Fourier Transform (IAAFT) method [47]. After that, measurements were compared to surrogates using one-sample single-sided t-tests, while adjusting for multiple comparisons through the False Discovery Rate (FDR) method [48]. The Granger causalities across the subset of glomeruli and somata displaying significant differences were then organized in two vectors. These were finally compared using a paired two-samples t-test to check for asymmetry in the causality strength between the two directions, which would confirm analysis validity. The underlying causality graphs are provided individually as Supplementary Materials. Representative spatiotemporal features. a) Example of the raw fluorescence signal from a focal plane across the antennal lobe with glomeruli at the center (numbered) and somata at the periphery (shaded). The superimposed regions-of-interest are color-coded by their response to isoamyl acetate, with magenta for activation, violet for absent, orange for inhibition. b) and c) Average and half standard deviation of glomerular and somata responses in the three response classes.

A. RESPONSE FEATURES SUMMARY
From the structural patterns of the raw fluorescence signals (Fig.5a), a total of 225 glomeruli and 537 somata were identified, yielding an ensemble of 1800 and 4296 corresponding time series, respectively. Out of all the glomerular responses, 28% were classified as activation and 17% as inhibition. Out of the responding somata responses, 26% were classified as activation and 13% as inhibition.
Visual inspection of the grand-average traces for the automatically-classified responses confirmed a well-evident difference between activation, inhibition, and absence of a response, for both the glomeruli (Fig. 5b) and the somata (Fig. 5c), which was overall in line with previous FIGURE 6. Spatiotemporal variability. a) Cumulative spatial distribution of the responses on the focal planes that were chosen to dissect the antennal lobe always at a similar depth (hue denotes overlaps, scale as in Fig. 5a). b) and c) Odor-related response temporal variability (bright line: average, shaded area: half standard deviation) recorded, respectively, from a representative soma and glomerulus.
observations. Inhibition responses were considerably shallower than activations, particularly for the glomeruli (≈ 2% vs. ≈ 5%). The averaged glomerular activations had welldefined, sharp edges alongside a variable overshoot. The somata activations, when considered individually, were also sharp; however, their timing was less consistent, leading to an overall smoother average [19], [49].
The spatial distribution of the responses did not reveal an anatomic segregation because most locations displayed both activation and inhibition as a function of odor type (Fig. 6a). Furthermore, considerable odor-related dynamical variability of the responses and associated classes could be observed within the individual glomeruli ( Fig. 6b) and soma (Fig. 6c). As the underlying relationship was unknown, for the purposes of this work, all time series were pooled. This implies that the transfer functions obtained represent the types of input-output transformations that the system can produce at a general level, as an ensemble average without reference to specific stimuli.

B. TRANSFER FUNCTIONS FROM ODORANT CONCENTRATIONS TO GLOMERULI
We first sought to determine the transfer functions between the external input signals to the olfactory system, namely the time-varying odorant concentrations, and the responses of its first processing stage, namely the glomeruli. For conciseness, in the below the response clusters are referred to with (x, o), where x = {A, I} denotes activation or inhibition and o = {P, Pt, Tp, T} denotes phasic, phasic-tonic, tonic-phasic or tonic output.
Responses within the activation and inhibition classes featured an inherent variability introduced by multiple factors, comprising different neural tunings across the glomeruli, odor-specific effects, and interactions between these aspects. As specified above, a further step of k-means clustering based on Euclidean distance, with k = 4, was performed to disentangle this heterogeneity at least partially. The resulting clusters for the activation responses revealed a continuum between an almost perfectly phasic, that is, derivative-like, response ( Fig. 7a) and an almost perfectly tonic, that is, proportional-like, response (Fig. 7d). This finding indicates that some components of the first processing stage react preferentially to the onset of an odor, whereas others encode its sustained presence. This diversification represents a highly conserved property of neural sensory architectures across species and modalities, serving, for example, to emphasize environmental changes [50], [51]. A similar observation could be established for the inhibition responses (Fig. 8), although their diversity was less evident, plausibly also owing to the smaller amplitudes.
The obtained transfer functions are reported in Table 1 for the activation and Table 2 for the inhibition responses. For the activation class, the majority of responses were determined to be intermediate between phasic and tonic, while the fit quality Q was appreciably higher for the tonic-phasic and tonic responses (Table 1). Differently, for the inhibition class, there was a preponderance of tonic-phasic responses, and the fit quality tended to be lower ( Table 2). The following features could be generally noted.
First, despite the inherent complexity of the system under investigation and the nonlinear, adaptive mechanisms it is based upon, these low-order linear models can reproduce the glomerular output signals elicited in response to the odorant presentation. Besides the fit quality measure Q, this could be visually ascertained by comparing the measured waveforms with the model outputs (Figs. 7 and 8, 1 st column). In all cases, the transfer function order was limited to 3, rendering these models particularly simple and computationally undemanding as compared to simulating a spiking neural network.
Second, the transfer function characteristics differ substantially between the clusters. Since, in the present context, the focus is on processing a stimulus having a welldefined onset, their features are conveniently illustrated with reference to canonical inputs, i.e., a Dirac (i.e., infinitelynarrow) impulse and a Heaviside unit step. For the activation responses (Fig. 7), all four identified systems are stable, such that the impulse response converges to zero, although  . Simulated model outputs in response to an impulse input (2 nd column) and a unit step input (3 rd column), alongside frequency response represented by the amplitude Bode plot (4 th column). All y-axes denote unitless amplitude, expressed in decibel for the Bode plots. Corresponding transfer functions given in Table 1. All step responses converge to different steady-state values, given that the models feature different static gains S 0 (Fig. 7, 3 rd column). Their diversity can be appreciated particularly by considering the relative overshoot, which, for perfectly tonic and phasic responses, by definition is σ = 0 and |σ | = ∞, respectively. Notably, the response for cluster (A,P) has a negative static gain S 0 < 0 and a remarkably narrow peak of positive response, leading to a . Simulated model outputs in response to an impulse input (2 nd column) and a unit step input (3 rd column), alongside frequency response represented by the amplitude Bode plot (4 th column). All y-axes denote unitless amplitude, expressed in decibel for the Bode plots. Corresponding transfer functions given in Table 2.
Considering next the frequency responses of these systems, clusters (A,P), (A,Pt) and (A,Tp), having a noticeable phasic contribution, are associated with band-pass behavior (Fig. 7ac, 4 th column), whereas cluster (A,T), having predominantly tonic output, is associated with low-pass behavior (Fig. 7d, 4 th column). Although the qualitative frequency response is the same across clusters (A,P), (A,Pt) and (A,Tp), the location and width of the passband are diversified, together with the level of attenuation at low frequencies. However, across all of them, the band between 1 and 10 rad/s appears similar. Altogether, these features assert that the input signal processing is qualitatively distinct between the clusters.
The transfer functions representing the inhibition responses delineate a comparable situation. All systems are stable, but their impulse and step responses display different characteristics (Fig. 8, 2 nd and 3 rd columns). Similar to the activations, all transfer functions have a negative and real pole alongside a pair of complex and conjugate poles, and the zeros are diversified ( Table 2). The most phasic response, for cluster (I,P), has two real and negative zeros, whereas the others, for clusters (I,Pt), (I,Tp) and (I,T), all have complex conjugate zeros.
The static gain and the damping factor of the poles vary between the clusters, such that the step response attains different steady-state values and the oscillations appear more persistent for cluster (I,Pt) than the others. Notably, the maximum overshoot is considerably lower than for the activations, namely σ = 1.75 vs. |σ | = 9.53, indicating a less marked differentiating action. Nevertheless, the gradient between clusters is clearly observable, in this case reaching σ = 0. The response time scales are similar, except for cluster (I,T), which features a longer rise time τ rise = 0.34 s.
Again, the frequency responses of these transfer functions are band-pass for clusters (I,P), (I,Pt) and (I,Tp) and lowpass for cluster (I,T), which reassures about the validity of the analyses. However, for cluster (I,T), due to the presence of a pair of weakly damped complex and conjugate zeros, a marked notch is apparently found around 10 rad/s; this should be interpreted cautiously.

C. INFORMATION FLOW BETWEEN GLOMERULI AND SOMATA
Over the time series pooled over all the bees and odors, the Granger causality analysis detected 19660 significant functional connections from the glomeruli to the somata, and 19447 connections in the opposite direction, with an underlying average graph density of ≈ 31% for 52±14 nodes. It should be again underlined that no univocal correspondence between glomeruli and somata could be established across individual bees, because of anatomical variability and the absence of means to standardize the identification of the latter. Therefore, the corresponding matrices could not be averaged together and are not shown here; the connectivity information was, in other words, only utilized within each set of time series to infer the probable underlying relationships between glomeruli and somata.
In line with the expectations for a primarily bottom-up, that is, feed-forward architecture, on average the causality strength was higher from the glomeruli to somata than in the opposite direction (two-samples t-test, t(39105) = 14.6, p < 0.001, 0.0188 ± 0.0154 vs. 0.0167 ± 0.0117, mean±standard deviation; Fig. 9a). The corresponding histograms overlapped substantially, and the distribution was right-skewed, more markedly so in the direction from the glomeruli to the somata. Namely, above a strength of G ≈ 0.025, the frequency of functional connections was consistently higher, for all bins, in the direction from the glomeruli to the somata than vice-versa (Fig. 9b,c).
This analysis could not fully discern the underlying connections, which are expected to be mainly feed-forward towards the somata [1]- [8]. Nevertheless, through comparing the measured to the surrogate time series at a statistical level, it was sufficiently capable of identifying the relevant combinations to support the mapping of transfer functions from the glomeruli to the somata. This was deemed acceptable for the present proof-of-concept study, focused on the transfer function-based modeling of these two processing stages. However, the elevated uncertainty associated with this analysis should be acknowledged.
The individual clusters identified at the glomerular level were selected but, to improve statistical reliability, only those glomerulus-soma pairs belonging to the same activation or inhibition class and having a Granger causality strength G > 0.025 were considered in estimating the transfer functions. At this threshold, there were 4429 causal interdependencies in the direction from the glomeruli to the somata, and 3643 in the opposite direction. Although limited to ≈ 20%, this asymmetry suggested that the inferred directed connections from the glomeruli to the somata could be utilized to tentatively reconstruct the transfer functions between them. At more stringent thresholds, the asymmetry increased drastically, with 793 vs. 338 connections (≈ 130%) for G > 0.05, and 86 vs. 9 connections (≈ 850%) for G > 0.1. While this confirms the validity of the analysis, in light of the fact that the number of connections dropped drastically, thus limiting statistical sampling, we opted for the more parsimonious threshold at G > 0.025. For the avoidance of doubt, only connections from the glomeruli to the somata were considered in the below, and these asymmetry indications are only reported as an indirect evidence of the putative validity of Granger causality results.

D. TRANSFER FUNCTIONS FROM GLOMERULI TO SOMATA
Having tentatively established the connectivity from the glomeruli to the somata, the transfer functions between them could be estimated. The classification into four clusters of the glomerular and somata responses gives rise to sixteen possible pairings of input/output response types, all of which could be observed in the recordings. The following analysis focused on the two extremes to reduce the amount of data and focus on the most representative pairings out of the four activation-related glomerular clusters. Namely, we considered the clusters having the most phasic (Fig. 7a) and the most tonic profiles (Fig. 7d), and, more specifically, we focused on the cases wherein they are paired with the most phasic and the most tonic profiles of somata responses. Four pairings were obtained from the data pooled in this manner, herein referred to with (x,i→o) where x = {A, I} denotes activation or inhibition and i, o = {P, T} denote phasic or tonic input and output. For each of them, a corresponding transfer function was derived, considering as input the average signal of a glomerular cluster and as output the average response of a somata cluster (magenta traces in Fig. 10, 1 st column). The same approach was applied to the inhibition responses (magenta traces in Fig. 11, 1 st  column). Because of the uncertainty associated with the indirect nature of the pairing between glomeruli and somata, as reflected in the relatively modest causality asymmetry, we parsimoniously refrain from presenting differences in prevalence between the response pairings.
The transfer functions obtained for the activation and inhibition responses (Tables 3, 4) confirm that third-order linear systems could also successfully capture the inputoutput behavior of the inferred glomeruli-somata couplings. This is evident considering the high overlap between the recorded somata outputs and the responses calculated using the corresponding estimated transfer functions (Fig. 10, 1 st column and 11, 1 st column). The fit quality Q of the activation-related transfer functions tended to be more consistently high as compared to the transformation of odorant concentrations into glomerular responses in the previous stage. Even more evidently in this instance, the inhibitionrelated transfer functions were associated with a markedly poorer fit quality, plausibly due to their lower amplitude.
A well-evident diversification of the responses is also appreciable in this case. However, considering firstly the activations (Table 3), it is noteworthy that the transfer functions have all real, and in particular negative, poles and zeros. Consequently, no oscillations appear in the impulse (Fig. 10, 2 nd column) and unit-step responses (Fig. 10, 3 rd  column).
As previously noted, the diversity of the transfer functions could be appreciated mainly based on the step responses, by considering the relative overshoot σ alongside the response time scale represented by the rise time τ rise . . Temporal curves (1 st column) for input (magenta), recorded output (gray), and model output (black). Simulated model outputs in response to an impulse input (2 nd column) and a unit step input (3 rd column), alongside frequency response represented by the amplitude Bode plot (4 th column). All y-axes denote unitless amplitude, expressed in decibel for the Bode plots. Corresponding transfer functions given in Table 3. . Temporal curves (1 st column) for input (magenta), recorded output (gray), and model output (black). Simulated model outputs in response to an impulse input (2 nd column) and a unit step input (3 rd column), alongside frequency response represented by the amplitude Bode plot (4 th column). All y-axes denote unitless amplitude, expressed in decibel for the Bode plots. Corresponding transfer functions given in Table 4.
For the pairings (A,P→P) and (A,T→T), the overshoot is shallow, with σ ≈ 0, and the response is fast, with τ rise ≈ 0.2 s, implying that the response is predominantly proportional, therefore propagating the input signal largely unaltered in its time course. On the other hand, for pairing (A,P→T), σ = 0 and the response is markedly slow, with τ rise ≈ 8.5 s, implying, over the relevant time scales, an integrative action that transforms a phasic input into a tonic output. Contrariwise, for pairing (A,T→P), the overshoot is substantial, with σ ≈ 40 and an almost null static gain S 0 ≈ 0, alongside a fast response τ rise < 0.1 s, implying a strongly differentiating action which transforms a tonic input into a phasic output.
For the frequency responses (Fig. 10, 4 th column), pairings (A,P→P) and (A,T→T), yielding respectively a phasic output in response to a phasic input and a tonic output in response to a tonic input, have similar characteristics, which resemble a low-pass filter featuring an inconspicuous resonance at ≈ 1 rad/s. The transfer function yielding a tonic output in response to a phasic input, pairing (A,P→T), also has low-pass filter behavior, however, without any peak. On the other hand, for pairing (A,T→P), yielding a phasic output in response to a tonic input, a strongly band-pass filter behavior is found, further confirming the strongly differentiating action at low frequencies.
Considering next the inhibition responses (Table 4), a different situation is found, which resembles the dissociation previously noted for the transfer functions from odorant concentrations to glomerular responses. Considering the markedly lower signal amplitude, these estimations and the associated phasic/tonic classifications should be considered as tentative (11, 1 st column), as also indicated by their relatively low quality Q. One transfer function, for pairing (I,T→T), has real and negative poles and zeros, whereas the others feature a pair of complex and conjugate poles with low damping alongside a pair of complex and conjugate zeros, TABLE 3. Transfer functions modeling the activation responses from the glomeruli to the somata. Q ∈ [0, 1] denotes fitting quality, S 0 and S max the static and maximum gains from the step response, σ the corresponding relative overshoot and τ rise the rise time. Corresponding plots shown in Fig. 10. all with negative real part. Consequently, damped oscillations could be observed in the impulse and unit-step responses (Fig. 11, 2 nd and 3 rd column). In particular, for the transfer function for pairing (I,T→P), yielding a phasic output in response to a tonic input, the damping factor is relatively small; this results in weakly-damped oscillations in the timedomain response alongside a marked peak in the frequency response (Fig. 11c, 3

rd and 4 th column).
Compared to the activations, in this case, the diversity as expressed by the relative overshoot σ and response time scale τ rise is less evident. In particular, all values of σ and τ rise are comparatively low, indicating that the integrative and differentiating action is weak. Nevertheless, the frequency responses of the pairings (I,P→P), (I,P→T) and (I,T→P) are characterized by low-pass filter behavior alongside a resonance peak having a different amplitude, reflecting the different damping in the system poles; on the other hand, for pairing (I,T→T), a band-pass filter behavior is observed.

IV. ANALOG SIGNAL PROCESSING A. CIRCUIT DESIGN
Next, we demonstrate how electronic circuits could be systematically designed based on these biologically-derived transfer functions. Importantly, circuit realization in this context should not merely be viewed as conducive towards hypothetical engineering applications. Instead, it should also, if not primarily, be intended as a way to confirm beyond the realm of numerical simulations the physical viability of the transfer functions obtained as models of mesoscopic neural dynamics.
For this purpose, we exemplify the calculation of the component values for the prototypical circuits shown in Fig. 12a and b. These realize the transfer functions corresponding, respectively, to all the pairings in Table 3 and to cluster (A,P) in Table 1. The arguments are generalizable to all other transfer functions considered in this work.
The circuits are based on active filters realized with operational amplifiers as the gain element, alongside resistors (R) and capacitors (C). They are obtained via cascading three (Fig. 12a) or five (Fig. 12b) stages. Each stage implements one of the multiplicative factors into which the corresponding transfer function can be decomposed, that is, a pole and a zero, or two complex and conjugate zeros together with two complex and conjugate poles. These circuits should be considered purely as examples, and other topologies could be equivalently selected for mapping a transfer function. They are, nevertheless, particularly parsimonious and therefore well-suited for discrete realization on a circuit board or miniaturized realization as part of a CMOS integrated circuit [52].
The first step in establishing a functional equivalence between the analog circuits and the transfer function models of the bee's olfactory response is to derive the transfer functions of the circuits shown in Fig. 12. This step can be completed by applying the standard methods of circuit theory [53]. For the circuit in Fig. 12a, one obtains the following transfer function, wherein R 6 and R 7 set the output gain (2) Since all the poles and zeros of this transfer function are real and negative, the same can be adopted to realize all the cases (A,P→P), (A,P→T), (A,T→P) and (A,T→T) in Table 3, once the component values have been determined according to the steps specified below.
On the other hand, the circuit in Fig. 12b implements the following transfer function The circuit comprises multiple stages, including a lowpass multiple-feedback configuration having a pair of complex and conjugate poles, two stages implementing a pair consisting of a negative pole and a negative zero, and two all-pass filters having a negative pole and a positive zero. It is necessary to select the poles and zeros such that i) the poles in −1/C 5 R 6 and in −1/C 8 R 9 cancel the negative zeros in −1/C 4 R 5 and in −1/C 7 R 8 , and ii) the pole in −1/C 6 R 7 is at a frequency sufficiently higher than the relevant band of the biological system so as to have a negligible impact. Then, the circuit can be used, for instance, as a physically-viable approximation of the transfer function for cluster (A,P) ( Table 1).  Fig. 12a). b) Third-order filter associated with a transfer function having a negative and real pole alongside a pair of complex and conjugate poles (circuit in Fig. 12b).
Having established the circuit topologies, the component parameters can be determined to realize the desired dynamics on the same time scale as a biological case. To illustrate this step, we return to pairing (A,P→P) in Table 3 and, for convenience, we recall the transfer function model then equate term-by-term Eq. (4) with Eq. (2), yielding Considering now the transfer function of cluster (A,P) in Table 1, that is, the one given by following the design guidelines described above, setting the high frequency pole at −1/ (C 6 R 7 ) = −k rad/s, and then adopting similar steps to the derivation used in the previous example, equating term-by-term Eq. (5) with Eq. (3), yields

B. CIRCUIT REALIZATION AND MEASUREMENT
Finally, we physically built and measured the circuits realizing the transfer functions corresponding to all the pairings in Table 3 ( Fig. 13a) and to cluster (A,P) in Table 1 VOLUME 10, 2022 TABLE 4. Transfer functions modeling the inhibition responses from the glomeruli to the somata. Q ∈ [0, 1] denotes fitting quality, S 0 and S max the static and maximum gains from the step response, σ the corresponding relative overshoot and τ rise the rise time. Corresponding plots shown in Fig. 11.. (Fig. 13b). For this purpose, off-the-shelf through-hole resistors (tolerance 5%) and electrolytic capacitors (tolerance 20%) were assembled onto a plug-board, together with operational amplifiers (type TL082, ST Microelectronics SpA, Agrate Brianza MI, Italy). The dual power supply voltage was set to ±9 V. Signal generation and digitization were performed through a USB-connected AD/DA converter (type AnalogDiscovery 2; Digilent Inc., Pullman WA). The waveforms input to the transfer functions (magenta traces in Fig. 7a and Fig. 10a-d) were upsampled by a factor of 100, then converted to analog signals at 1 kSa/s, 14-bit. To prevent saturation of the intermediate stages, the peak amplitude was set to about 100 mV. The signals generated by the circuits were digitized at 800 Sa/s, 14-bit, downsampled, and re-aligned to the biologically-recorded outputs (orange trace in Fig. 7a and gray traces in Fig. 10a-d). To assess the degree of similarity, the point-to-point linear correlation was calculated. It is worth underlining that, given the fact that the signal amplitude swings at all stages were small compared to the power supply voltage and that the input signals were very slow compared to the gain bandwidth of the chosen integrated circuits, the operational amplifiers could be treated as ideal entities. In other words, insofar as these assumptions are met, according to fundamental circuit theory one can abstract from the non-linear and non-ideal characteristics of the underlying transistors and simply treat each stage as a purely linear entity [52]. The underlying recordings are provided as Supplementary Materials.
Because off-the-shelf resistors and capacitors values are discretized according to standard tables, the component choices were adjusted to approximate the intended parameter settings as closely as possible. To minimize drifting which would otherwise cause saturation, a finite static gain needs to be set. For that purpose, 470 k resistors were added in parallel to the RC feedback networks, namely C 3 /R 3 and C 5 /R 5 in Fig. 12a, and C 4 /R 5 and C 7 /R 8 in Fig. 12b. Together with resistor tolerances, this introduced minor mismatches in the overall gain. For the circuit in Fig. 12b, to avoid an excessively high gain in the fourth stage, we set k = 20; while this value is comparatively low, the corresponding high-frequency pole had a negligible impact. Moreover, to maximize the signal-to-noise ratio during measurements, the input signal was amplified ×10, and the output was accordingly divided (not shown). At this point, it is opportune to point out a difference with respect to the initial scenario. In the biological recordings, the signal-to-noise ratio intended as the random dispersion of the evoked responses was on the order of unity, as visible in Fig. 5b. In the realized circuits, the same was approximately 100-fold higher. In other words, the noise present in the initial recordings was treated merely as a source of parametric uncertainly in the transfer functions, and not replicated. While this is unavoidable in this linear analysis, it should be borne in mind that the initial noise is not only due to instrumental and physiological contamination, but stochastic dynamics form an essential aspect of neural behavior also at the mesoscopic scale [1], [3], [10], [13]. Such dynamics would not in any case have been realistically reproduced by simply choosing a lower signal-tonoise arrangement of the linear electronic circuits, therefore, we parsimoniously aimed to maximize the signal-to-noise ratio while explicitly acknowledging this simplification.
The circuit-generated output, visible in Fig. 14b, also had a near-complete overlap with the corresponding biological recording (r = 0.99, p < 0.001).
The circuit-generated output, visible in Fig. 14c, also had a near-complete overlap with the corresponding biological recording (r = 0.99, p < 0.001).
The circuit-generated output, visible in Fig. 14d, also had a near-complete overlap with the corresponding biological recording (r = 0.97, p < 0.001).
For cluster (A,P) in Table 1 .
Also in this case, the circuit-generated output, visible in Fig. 14e, had near-complete overlap with the corresponding biological recording (r = 0.99, p < 0.001). The five realized circuits, therefore, reproduced with remarkable accuracy the neural responses, despite the approximations associated with the finite transfer function order and realization using discrete electronic components.

A. TRANSFER FUNCTIONS AS A ''GRAY BOX'' MODEL
This study was predicated on the ability to simultaneously record the neural responses during olfactory stimulation in the glomeruli and somata using two-photon calcium imaging. This enabled determining models of biological signal processing across the first two stages of the honey bee's olfactory system. The first processing step in the antennal lobe can be dissected into two input-output subsystems, namely, from the ORN dendrites to a glomerulus and from the glomerulus to the PN somata. Therein, channels are laterally interlinked by LNs, which were found in the fruit fly to be both excitatory and inhibitory [15], [16]. In the honey bee,  Table 3 (Eq. 9), e) cluster (A,P) in Table 1 (Eq. 10).
thus far, only inhibitory local neurons have been identified directly [54].
According to established models of the antennal lobe, one would expect limited changes in the signal propagating from the PN dendrites in the glomeruli to the somata, since the coupling between PNs via local neuron is assumed to be limited to the glomerular region [11]. However, calcium imaging paired with electrophysiological recordings in the moth AL suggests that the PN dendritic calcium signal from the glomeruli tracks the postsynaptic potential, whereas the calcium signal in the somata is a good measure of the PN spiking output [55]. Moreover, the PN dendritic signal in the glomeruli seems to depend primarily on acetylcholine channel influx rather than intracellular potential, whereas the somata signal depends more on voltage than acetylcholine [56]. Despite these intricacies, the inferences about olfactory coding that can be drawn from calcium imaging and direct electrophysiological recordings are essentially overlapping [57].
The sheer difficulty in ascribing the diversified responses across the glomeruli and somata to well-defined architectural VOLUME 10, 2022 and physiological aspects ultimately motivated us to consider an ''agnostic'' approach. Namely, this entailed capturing at the mesoscopic scale the changes in the input signals across these stages, without developing a detailed model of the system internal behavior, accepting instead a socalled ''gray box'' representation in terms of a low-order transfer function [24], [53]. On the one hand, this approach is unrewarding in terms of attaining a detailed understanding. However, on the other, it is very effective in yielding efficient and compact mathematical or even physical artificial replicas of biological processes. These can be used for a multitude of purposes, as exemplified by the active filter circuits demonstrated in this work [31], [32]. Such a situation is representative of the role of ''gray box'' models across diverse areas of engineering and science [25], [26].

B. GLOMERULAR RESPONSES
A principal finding from the outputs at the glomerular level is that the activity of the PNs shows a broad range of excitation and inhibition responses, with widely varying dynamics, particularly regarding transiency. A classification based on response parameters did not clearly correspond to glomerulus location and evoking odors, suggesting an intricate interaction among manifold mechanisms of signal transduction. In general, the PN response profiles were found to exhibit the typical phasic-tonic shape of activation responses, showing a steep rise and a relatively rapid accommodation (Fig.7) after stimulus end; furthermore, a phasic undershoot followed most responses. Conversely, inhibition responses were, in the majority of cases, found to be tonic (Fig.8). These response types have been consistently observed in previous studies of PN activity in the honey bee [58], and can be accounted for by lateral inhibition via the LNs setting in shortly after the initial excitation.
It has been posited that these transformations broaden the PN tuning, producing more uniform distances between odor representations in the PN coding space [59]. Here, cluster analysis of the responses showed that the diversifying trait is the different weighting between phasic and tonic contributions, ranging from purely spike-like to almost box-like responses. For convenience, this was illustrated using a four-level categorization, though the results are not textcolorredcritically dependent on this choice. This diversification is plausibly introduced by coupling to the inhibitory network of local neurons, whose contribution varies in number and strength. Spike-like responses foremost emphasize the odor onset times, whereas box-like responses reproduce the stimulus pulse shape and intensity [60].
In turn, this implies that the individual glomerular responses contribute differently to different coding mechanisms [61]. Rate or amplitude coding use information about the neuronal firing rate or equivalently the calcium amplitude to represent information about a stimulus [62]. This requires integrating glomerular responses over time and hereby delivers precise information on the stimulus intensity [63]. Conversely, the mechanism of latency or first-spike coding is based on detecting the onsets of single glomerular responses, wherein their relative time differences provide rapidlyavailable information about odor identity [22], [64].

C. SOMATA RESPONSES
Thanks to instrumental advances, in this study we could, for the first time, acquire the calcium signals in parallel from a large number of glomeruli-somata pairs. The results regarding their association need to be interpreted carefully since no structural connectivity information was available. Connectivity was, therefore, inferred indirectly via statistical estimates of Granger causality based on the recorded time series. Nevertheless, some tentative conclusions could be drawn. We confirm that responses in the somata are generally more transient than those in the glomeruli [19], [49]. In about half of the activation responses, phasic and tonic profiles retain their shapes (Fig.10 a,d). However, in the other half, a significant change in the response profile could be detected, such that the majority of output signals are more spikelike compared to the corresponding inputs (Fig.10c vs. b).
Remarkably, for the inhibition responses, the tendency was the opposite; although, also in this case, about half of the output signals did not significantly change shape (11 a,d), those that did show a preponderance of prolonged responses (11 b with respect to c).
Indirectly, this situation suggests that, in a fraction of cases, a further inhibitory lateral modulation may be acting on the signal post-synaptically along its way from the PN dendrites to the neuronal body/axonal cone. While, to the authors' knowledge, this interaction is not included in the standard antennal lobe models, it has been proposed before, predicated on experiments showing pulse shape transformations similar to the ones observed here. For example, a reduction in the excitation pulse length at the level at the somata was reported in the moth [55]. The reason for such a transformation is plausibly to be found in the transition between coding schemes, namely, from dense coding in the antennal lobe to sparse coding in the mushroom body, which unavoidably also requires a sparsification of pulses in time [65]. Consistent results for inhibition responses were found in the moth [66]. Whether the occasionally-observed opposite transformation, namely prolonged excitatory and shortened inhibitory responses, could also be explained by lateral coupling or not, depends on the presence of local excitatory neurons. These have been found in the fruit fly [15], [16] but, thus far, not in the honey bee. However, the broad tuning of PN responses in the honey bee AL is an argument for their existence [60].
An alternative and equally plausible explanation is that these transformations are not describing processes in single neurons, but that the inputs and outputs come from highly correlated yet distinct cells. Clearly, it also cannot be excluded that these transformations may arise as an artifact of the selection of glomeruli-somata pairs based on Granger causality. That appears plausible in light of the limited asymmetry observed. Future work should address this aspect, for example, as a function of different causality thresholds to sparsify the graphs. Nevertheless, the wide variety of output signal shapes found in our recordings appears well in line with calcium imaging studies performed downstream, at the axon terminals (buttons) of the mushroom body calyces; at that point, however, further modulation by feedback neurons from the mushroom body plausibly also plays a role [67].

D. PUTATIVE TRANSFORMATION MECHANISMS
The observed relationships between the glomerular and somata responses, as said, point to the possible action of postsynaptic inhibition beyond the glomeruli.This would realize a form of adaptive gain control [55], which may underpin the observation that second-order neurons display broader tuning and more complex responses [14]. Namely, PN responses rise and accommodate rapidly, emphasizing odor onset, such that weak ORN inputs are amplified, whereas strong inputs are not. This nonlinear transformation broadens PN tuning, yielding more uniform distances between odor representations in the PN coding space [59]. However, these gain effects can not be directly inferred from the measured amplitudes of the calcium-induced fluorescence changes (Fig. 5,6). As stated above, on the one hand, the response amplitudes depend on uncontrollable factors such as dye uptake and, on the other, although the measurements come from the same cell, the biophysical contributions to the calcium concentrations vary strongly between dendrites and soma [55], [56]. This resembles the situation for other indirect measures of neural activity whose coupling functions are underdetermined, such as the blood oxygen signal in functional MRI and near-infrared spectroscopy [68]. Therefore, electrophysiological measurements are necessary to make quantitative assertions about gain and therefore fully determine neuronal transfer functions, as extensively exemplified and argued in Refs. [27], [42]. The present results inform how to target such measurements in future work.
From another perspective, the transformation from ORN to PN described by the normalization model equalizes magnitudes and decorrelates responses to different odors through feedforward nonlinearities and lateral suppression. This allows highly selective and sparse responses to be generated through the application of a fixed threshold, for which response magnitude equalization is crucial [69]. This processing step may implement a spatiotemporal sparsification through spike frequency adaptation [65], and more generally a central role has been posted for self-organization and criticality in olfactory processing [2]. Altogether, these arguments unavoidably raise the question of whether modeling this knowingly complex processing via transfer functions has any true explanatory value.
There are, in particular, three issues to be considered. First, the transfer functions are linear, which is at odds with manifold evidence of nonlinear processing. Second, they embody a pure feed-forward view of the processing architecture, which is also at odds with the view that lateral inhibition plays a key role. Third, their parameters are assumed to be fixed, whereas, in biological systems, dynamical activity and adaptiveness imply continuous and pervasive retuning.

E. THE SCOPE FOR TRANSFER FUNCTIONS
The relevance or otherwise of the proposed approach based on determining transfer functions from neurophysiological time series, then synthesizing analog electronic circuits based on them, critically depends on the purpose of investigation. Evidently, this approach cannot replace the microscale, bottom-up analysis of neural microcircuits, which is the basis of the field of spiking neuromorphic circuits. That approach yields analog CMOS architectures which mimic individual neurons and their interconnections, while having circuit variables explicitly associated to physiologicallymeaningful quantities, such as ion currents [70]. Similarly, this approach cannot replace large-scale computational investigation wherein complex dynamical phenomena can be flexibly and extensively explored, aiming, for example, to explain self-organization phenomena and other emergent properties [2]. Plausibly, transfer functions would also be incapable of accounting for the myriad of intricate interrelationships between odor type and spatial distribution of neural activity, which have been outright neglected in this work. As considered in the previous Section, by their nature, they represent a highly simplified, linear and deterministic account of neural dynamics, which by construction neglects both their stochasticity and their non-linear nature. Consequently, and no less importantly, transfer functions as used in this work represent a ''snapshot'' of the inputoutput relationship. Since the coefficients cannot change as a function of the previous inputs, they are inherently unsuitable for creating hardware systems capable of learning. As such, they are just models of signal filters, rather than adaptive devices.
Instead, the purpose here was to show that a minimalistic, high-level empirical representation of the mesoscale-level system responses can be agnostically obtained using transfer functions, that is, choosing ab initio a ''gray box'' approach which prioritizes compactness and ease of realization. For example, in a biological scenario, the conversion between phasic and tonic responses may involve self-tuning microcircuits considerably more complex than these tentative linear approximations [71], [72]. Nevertheless, in some cases a detailed representation of such mechanisms may not be necessary to reproduce their overall response satisfactorily. In particular, this level of complexity reduction may be the key to unlocking the practical feasibility of bio-inspired engineering in real-world applications.y The empirical observation of similar responses should clearly not be misinterpreted as implying a structural correspondence between the electrical components in the circuits (e.g., resistors and capacitors) and elements of neural circuitry. Such a one-to-one correspondence cannot be established by the present kind of research. In fact, it cannot be established at all based on two-photon microscopy and necessitates methods such as single-cell recordings. These models should, therefore, be intended as phenomenological rather than structural, and as mesoscopic rather than microscopic representations.
At the same time, an important limitation of the present work which should be acknowledged is the fact that no attempt was made to decode the complex association between the type of odor presented and the spatiotemporal distribution of activation over the glomeruli and somata [6]- [10]. Pooling all the data and resorting to a fundamental and linear tool such as Granger causality to establish the correspondence between the glomeruli and the somata may be viewed as an abdication from answering this question. However, the focus of this paper was on delivering a compact and linear account of signal processing [24], [44], [45]. On the one hand, such a step was necessary because, while the anatomical variability of the glomeruli is limited, that of the somata is considerably larger due to their sheer number and small size [11], [21], [22], [34]. On the other hand, we acknowledge that future work should address the application of nonlinear tools from machine learning, such as deep neural networks, to decoding this representation, eventually informing a more detailed understanding of the entire information processing chain.
Despite all the conceptual limitations mentioned above, the accuracy of the synthesized analog circuits in reproducing the initial recordings was striking, with a near-perfect correlation observed between the biological and electronic experiments. From this, one concludes that, at the mesoscopic scale, regardless of the underlying complexity, even a low-order linear approximation can be adequate to account for the ensemble responses around a particular operating point. On the one hand, the proposed transfer functions are unlikely to be able to account for concentration-related and odorspecific effects, which were not addressed in this work. On the other, the obtained circuits are incomparably simpler to realize than complex neuromorphic architectures, and their parameters can be determined straightforwardly via closedform equations. Equally, by their nature, they are orders of magnitude less power-consuming than the numerical simulation of any realistic model. As regards neuromorphic engineering in particular, it is worth remembering that, even though presently the focus is mainly on single neuronlevel replication, in the past architectures based on capturing at the mesoscale the behavior of entire populations of neurons, such as slow fluctuations in local field potential in response to light, have been proposed and successfully deployed [73]. The scope of the present transfer functionbased modeling should therefore be considered in that regard, and particularly with an eye to the practical realization of built-in sensor front-ends, such as for artificial noses. In this sense, the ability to synthesize a biologically-plausibly response with a closed-form approach, based on standard operational amplifiers, is a notable advantage over the known intricacies of iterative design and parameter tuning associated with spiking neuromorphic circuits [31], [32].
[71] C. G. Galizia   He has authored more than 150 articles and five patents. His research interests include non-linear dynamical systems, chaotic oscillators, reconfigurable analog and digital computing, functional magnetic resonance imaging, advanced techniques for bio-signal analysis, brain-machine/computer interfaces, and robotics. He is also a member of the Institute of Electronics, Information, and Communication Engineers (IEICE), the Institute of Electrical Engineers of Japan, and the Japan Neuroscience Society. He is also Chartered Engineer, a European Engineer, and a member of the Institution of Engineering and Technology, U.K, a Chartered Physicist and a member of the Institute of Physics, U.K., and a Chartered Scientist and a member of the Institute of Physics and Engineering in Medicine, U.K. SIMONE MONACHINO received the bachelor's degree in biological sciences from the University of Rome Tor Vergata, Italy, in 2017, the master's degree in computational neuroscience and neuroinformatics from Newcastle University, U.K., in 2018, and the master's degree in cognitive science from the University of Trento, Italy, in 2020. His research interests include computational neuroscience, neuroinformatics, and the mechanisms of neural coding. He is a member of the Organization for Computational Neurosciences (OCNS), the Interna-