A Molecular Odorant Transduction Model and the Complexity of Combinatorial Encoding in the Drosophila Antenna

In the past two decades, a substantial amount of work characterized the odorant receptors, neuroanatomy and odorant response properties of the early olfactory system of Drosophila melanogaster. Yet many odorant receptors remain only partially characterized and, the odorant transduction process and the axon hillock spiking mechanism of the olfactory sensory neurons (OSNs) have yet to be fully determined. Identity and concentration, two key aspects of olfactory coding, originate in the odorant transduction process. Detailed molecular models of the odorant transduction process are, however, scarce for fruit flies. To address these challenges we advance a comprehensive model of fruit fly OSNs as a cascade consisting of an odorant transduction process (OTP) and a biophysical spike generator (BSG). We model identity and concentration in OTP using an odorant-receptor binding rate tensor, modulated by the odorant concentration profile, and an odorant-receptor dissociation rate tensor, and quantitatively describe the ligand binding/dissociation process. We model the BSG as a Connor-Stevens point neuron. The resulting combinatorial encoding model of the Drosophila antenna provides a theoretical foundation for understanding the neural code of both odorant identity and odorant concentration and advances the state-of-the-art in a number of ways. First, it quantifies on the molecular level the combinatorial complexity of the transformation taking place in the antennae. The concentration-dependent combinatorial code determines the complexity of the input space driving olfactory processing in the downstream neuropils, such as odorant recognition and olfactory associative learning. Second, the model is biologically validated using multiple electrophysiology recordings. Third, the model demonstrates that the currently available data for odorant-receptor responses only enable the estimation of the affinity of the odorant-receptor pairs. The odorant-dissociation rate is only available for a few odorant-receptor pairs. Finally, our model calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies.


Introduction
The odorant response of olfactory sensory neurons (OSNs) in the Drosophila antennae has been experimentally characterized by multiple research groups (de Bruyne et al., 2001;Hallem and Carlson, 2006;Smart et al., 2008), and their results combined into a single consensus database, called the DoOR database (Galizia et al., 2010;Münch and Galizia, 2016). A key functionality of OSNs is to jointly encode both odorant identity and odorant concentration (Masek and Heisenberg, 2008;Yarali et al., 2009;Nagel and Wilson, 2011;Wilson, 2013;Martelli et al., 2013). A single odorant stimulus usually activates multiple OSN groups expressing the same receptor type, while different odorants activate different OSN groups (Laurent, 1999;Hopfield, 1999;Touhara and Vosshall, 2009). The identity of an odorant is combinatorially encoded by the set of responding OSN groups (Malnic et al., 1999), and the size of OSN set varies as the concentration changes (Hallem and Carlson, 2006). The temporal response of an OSN simultaneously represents the information of odorant concentration and concentration gradient also known as 2D odorant encoding (Kim et al., 2011b(Kim et al., , 2015. These two aspects of olfactory coding, identity and concentration, originate in the odorant transduction process (Hallem et al., 2004;Nagel and Wilson, 2011). However, detailed molecular models of the odorant transduction process are scarce for fruit flies.
To address these challenges we advance a comprehensive model of fruit fly OSNs as a cascade consisting of an odorant transduction process (OTP) and a biophysical spike generator (BSG). We model identity and concentration in OTP by an odorant-receptor binding rate tensor, modulated by the odorant concentration profile, and an odorant-receptor dissociation rate tensor, and quantitatively describe the ligand binding/dissociation process. OSNs are distributed across the surface of maxillary palp and the third segment of antenna (de Bruyne et al., 2001;Couto et al., 2005;Song et al., 2012). Since there is no commonly accepted terminology in the literature for naming these two olfactory appendages as a single entity, and in order to avoid potential confusion, we will refer to the set of all OSNs on one side of the fly brain as an antenna/maxillary palp (AMP) local processing unit (LPU) (Chiang et al., 2011).
To biologically validate our modeling approach, we first propose an algorithm for estimating the affinity and the dissociation rate of an odorant-receptor pair. We then apply the algorithm to electrophysiology recordings and estimate the affinity and dissociation rate for three odorant-receptor pairs. Second, we evaluate the temporal response of the OSN model with a multitude of stimuli for all three odorant-receptor pairs. The output of the model closely reproduces the temporal responses of OSNs obtained from in vivo electrophysiology recordings (Kim et al., 2011b(Kim et al., , 2015 for all three odorant-receptor pairs across all three types of stimuli. Lastly, we evaluate the model at the OSN antennae population level. We first empirically estimate the odorant-receptor affinity using the spike count records in the DoOR database for 24 receptor types in response to 110 odorants (Münch and Galizia, 2016). With estimated affinity values, we simulate the temporal response of the OSN population to staircase odorant waveforms. The output of simulated OSN population demonstrates that the odorant identity is encoded in the set of odorant-activated OSN groups expressing the same receptor type, and, more importantly, the size of the set expands or reduces as the odorant concentration increases or decreases.
The fruit fly OSN model presented here provides a theoretical foundation for understanding the neural code of both odorant identity and odorant concentration. It advances the state-of-the-art in a number of ways. First, it models on the molecular level the combinatorial complexity of the transformation taking place in Drosophila antennae OSNs. The resulting concentration-dependent combinatorial code determines the complexity of the input space driving olfactory processing in the downstream neuropils, such as odorant recognition and olfactory associative learning. Second, the model is biologically validated using multiple electrophysiology recordings. Third, the OSN model demonstrates that the currently available data for odorant-receptor responses only enables the estimation of the affinity of the odorant-receptor pairs. Finally, our model calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies.

Model
Detailed biophysical models for the odorant transduction process have been proposed for worms and vertebrates. Such models are scarce for insects and, in particular, for fruit flies. Dougherty et al. proposed a frog odorant receptor model that exhibits a complex temporal response (Dougherty et al., 2005). Rospars et al. proposed a model that characterizes the steady state response of OSNs for rats and moths (Rospars et al., 1996(Rospars et al., , 2008. The model stands out for its simplicity and modeling clarity, while lacking temporal variability. Other notable models appeared in (Suzuki et al., 2002;De Palo et al., 2012). Recently, Cao et al. published a phenomenological model to characterize the peak and the steady response of sensory adaption for fruit fly OSNs (Cao et al., 2016). Gorur-Shandilya et al. proposed a two-state model for the fruit fly odorant receptors that can reproduce Weber-Fechner's law observed in physiological recordings (Gorur-Shandilya et al., 2017). In addition, De Palo et al. (De Palo et al., 2013) proposed an abstract/phenomenological model with feedback mechanism that characterizes the common dynamical features in both visual and olfactory sensory transduction. Except for the transduction current recorded for studying sensory adaptation (Cao et al., 2016), reproducing the temporal response of the AMP LPU on either side of the brain has been scarcely investigated in the literature. In particular, 2D odorant encoding has not yet been successfully modeled. To address these challenges, we model the OSNs as a cascade consisting of an odorant transduction process (OTP) and a biophysical spike generator (BSG), as shown in Figure 1. The OTP model consists of an active receptor model and a co-receptor channel model (Larsson et al., 2004;Benton et al., 2006). The BSG model we employ here is based on the Connor-Stevens point neuron model (Connor and Stevens, 1971). The spike trains generated by the BSGs contain the odorant identity, odorant concentration, and concentration gradient information that the fly brain uses to make odorant valence decisions. Figure 2. Three dimensional odorant-receptor binding rate tensor . For a given neuron = 1, 2, ..., , the binding rate values are denoted by [ ] , for all = 1, 2, ..., , and = 1, 2, ..., . For the fruit fly, the total number of neurons expressing the same receptor type is about = 25, and the total number of receptor types is around = 60. is the number of all odorants that the fruit fly senses.

The Odorant Transduction Process Model
Two research groups have published widely different results on the OR transduction process in fruit flies (Sato et al., 2008;Wicher et al., 2008;Kaupp, 2010). As the exact signaling of the transduction cascade in fruit flies is still inconclusive, our approach focusses here on constructing a minimal transduction model. Called the fruit fly odorant transduction process (OTP) model, it extends the model proposed by (Rospars et al., 1996(Rospars et al., , 2008 by incorporating the essential features of temporal dynamics of other computational models, such as the one proposed by (Dougherty et al., 2005), while at the same time exhibiting the calcium dynamics of (Cao et al., 2016). In the latter work, the temporal dynamics of fly's OSN vanish in the absence of extracellular calcium. Notably, the calcium dynamics considered here constitutes a feedback mechanism that is similar to but also different from the one in the abstract model proposed by (De Palo et al., 2013).
Olfactory transduction in fruit flies from airborne molecules to transduction current involves a number of steps (Kaissling, 2013;Leal, 2013): i) absorption of odorant molecules through the sensillum surface, binding between odorant molecules and odorant-binding proteins (OBPs), and diffusion of bound OBPs through the aqueous sensillar lymph to OSN dendrites, ii) odorant-receptor binding/dissociation, and iii) opening of ion channels that results in transduction current. The first step is known as the "peri-receptor" processing, the second step is referred as the bound receptor generator and the third step as the co-receptor channel. Taken together, they represent the fruit fly odorant transduction process.
We propose an olfactory transduction process model that consists of an active receptor model and a co-receptor channel model (Larsson et al., 2004;Benton et al., 2006). The active receptor contains a peri-receptor model and a bound-receptor model. The peri-receptor process is modeled as a linear filter that describes the transformation of an odorant concentration waveform as odorant molecules diffuse through sensilla walls towards the OSN dendrites. The bound-receptor model encodes odorant identity and odorant concentration with a binding rate tensor, modulated by the odorant concentration profile, and a dissociation rate tensor. The odorant concentration profile is defined as the linearly weighted sum of the filtered odorant concentration and the filtered concentration gradient. Modulation is modeled here as a product. The co-receptor channel represents the ion channel gated by the atypical co-receptor (CR), Or83b. The calcium channel models the calcium dynamics, and provides a feedback mechanism to the co-receptor channel. The transduction current generator models the transmembrane current through the co-receptor channel.

The Active Receptor Model
The fruit fly active receptor model quantifies the binding and the dissociation process between odorant molecules and odorant receptors. As introduced here, the model centers on the rate of change of the ratio of free receptors versus the total number of receptors [ 0 ] expressed by neuron : where [ 1 ] is the ratio of ligand-bound receptors. Here = 1, 2, ..., , is the receptor type, = 1, 2..., , denotes the odorant and = 1, 2, ..., , denotes the neuron index. In Equation 1 above the ratios [ 0 ] and [ 1 ] are entries of the 3D tensors 0 and 1 , respectively. The 3D tensor with entries [ ] is called the odorant-receptor binding rate and models the association if the RHS is positive and zero otherwise. The RHS is the weighted sum of the filtered odorant concentration and the filtered concentration gradient ∕ with [ ] denoting a weighting factor. The impulse response of the linear filter ℎ = ℎ( ), ∈ ℝ (ℝ denotes the set of real numbers), models the "peri-receptor" process that describes the transformation of odorant concentration waveform as odorant molecules diffuse through sensilla walls towards the dendrites of OSN (Fleischer et al., 2017). For simplicity, the dependence of ℎ( ) on the geometry of the sensillum and the diffusion of odorant molecules across the sensillar lymph is not considered here. ℎ( ) in Equation 2 is the impulse response of a low-pass linear filter that is usually defined in the literature in frequency domain. Alternatively, ℎ( ) is the solution to the second-order differential equation, see Appendix 1.
Note that the odorant transduction models in the literature only consider the odorant concentration but not the odorant concentration profile as the input to the transduction cascade (De Palo et al., 2012Rospars et al., 1996Rospars et al., , 2008Cao et al., 2016). As we will show, the odorant concentration profile is critical for modeling 2D odorant encoding.
We assume that receptors only have two states, either being "free" or "bound", i.e., [ 0 ] + [ 1 ] = 1. Then, Equation 1 amounts to, Equation 3 maps the input given by the product between the binding rate and the odorant concentration profile, and the dissociation rate, i.e., . In what follows this map will be called the bound-receptor generator.

The Co-Receptor Channel Model
The fruit fly co-receptor channel and a calcium channel appear in a feedback configuration. Each of these components has its specific functionality. The co-receptor channel represents the ion channel gated by the atypical co-receptor (CR), Or83b. The calcium channel models the calcium dynamics, and provides a feedback mechanism to the co-receptor channel.
Next, we walk through each of the three equations of the co-receptor channel model. The key variables involved in the proposed odor transduction model are summarized in Table 1.
(1) The rate of change of the gating variable of the co-receptor channel [ 2 ] : where 2 and 2 are scalars indicating the rate of activation and deactivation of the gating variable, respectively, and where 3 and 3 are scalars indicating the rate of increase and decrease of the state variable.
(3) Finally, the transduction current [ ] is given by: where and are scalars, and denotes the maximal amplitude of the current through the co-receptor channel, whose value is empirically determined through parameter sweeping.

Biophysical Spike Generator Model
We restrict our choice of the spiking mechanism of OSNs to biophysical spike generators (BSG) such as the Hodgkin-Huxley, the Morris-Lecar, and the Connor-Stevens point neuron models. For simplicity of presentation, we only consider in Appendix 2 the Connor-Stevens (CS) point neuron model.

Biological Validation of the OSN Model
The essential functionality of OSNs is to jointly encode both odorant identity and odorant concentration. To address these two functional aspects we modeled each OSN as an OTP/BSG cascade. To validate our approach, we examine here the response of the OSN model to odorant waveforms that were previously used in experiments with different odorants and receptors, and compare the model responses with electrophysiological recordings. We first estimate the affinity and dissociation rate for the (acetone, Or59b) pair using two different datasets of electrophysiology recordings. Second, we evaluate the temporal response of the Or59b OSN model to acetone with a multitude of stimuli, including step, ramp and parabola waveforms. We further interrogate the model with staircase and white noise waveforms. Our results show that the model closely matches the complex temporal response of Or59b OSNs from electrophysiological recordings. Lastly, we evaluate the affinity and dissociation rate for different odorant-receptor pairs including (methyl butyrate, Or59b) and (butyraldehyde, Or7a). Estimating the Affinity, Binding and Dissociation Rates for (Acetone, Or59b) We applied Algorithm 1 (see Methods and Materials section) to estimate the affinity, the dissociation rate, and the binding rate for (acetone, Or59b) by using two different datasets of electrophysiology recordings. The source of the two datasets is given in Appendix 3. Each of the two datasets contains the PSTHs obtained from the response of OSNs expressing Or59b to acetone step waveforms with different concentration amplitudes. As required by Algorithm 1, we first retrieved the peak and steady state spike rates from the PSTH in response to each concentration amplitude recorded in the datasets. Second, for each of the two datasets, we used the steady state spike rate to estimate the value of the affinity for each concentration amplitude, and computed the mean and variance of the affinity as shown in Figure 3.A. With the mean of the estimated affinity, we then used the peak spike rate to estimate the value of the dissociation rate for each concentration amplitude, and computed the mean and variance of the dissociation rate as shown in Figure 3.B. For the first dataset, the mean and variance of the estimated affinity are 3.141 ⋅ 10 −4 and (1.312 ⋅ 10 −4 ) 2 , respectively, and the mean and variance of the estimated dissociation rate are 1.205 ⋅ 10 1 and (3.900 ⋅ 10 1 ) 2 , respectively. For the second dataset, the mean and variance of the estimated affinity are 3.201 ⋅ 10 −4 and (1.001 ⋅ 10 −4 ) 2 , respectively, and the mean and variance of the estimated dissociation rate are 1.389 ⋅ 10 1 and (1.262 ⋅ 10 1 ) 2 , respectively. The values of the affinity estimated from the two datasets are almost identical, while the two estimated dissociation rates are marginally different. This is because the steady state spike rates of the two datasets are similar, but the peak spike rates of the two datasets differ slightly, as shown in Figure 1 in Appendix 3.
Evaluating the Temporal Response of the Or59b OSN Model to Acetone To evaluate temporal response, we stimulated the OSN model with multiple odorant stimuli that were previously used in experiments designed for characterizing the response to acetone of OSNs expressing Or59b. For all OTP models considered below we set the odorant-receptor binding rate to 3.141 ⋅ 10 −4 and the odorant-receptor dissociation rate to 1.205 ⋅ 10 1 .
Response of the Or59b OSN Model to Step, Ramp and Parabola Acetone Waveforms  (Kim et al., 2015) using the original raw data).
We first evaluated the response of the Or59b OSN model to step, ramp, and parabola stimulus waveforms as shown in the first row of Figure 4. The temporal response of the OTP/BSG cascade (the fourth row of Figure 4) is similar to the one of the OTP model (the third row of Figure 4). For step stimuli, the OTP/BSG cascade generates a chair-shaped response by first picking up the gradient of the concentration right after the onset of the odorant, and then gradually dropping down to a constant value, that encodes the step value of the amplitude. For ramp stimuli, the initial response of the OTP/BSG cascade rapidly increases, and then it plateaus as the gradient of ramp stimuli becomes constant. Lastly, for the parabola stimuli, the response of the OTP/BSG cascades resembles a ramp function, that corresponds to the gradient of parabola stimulus waveforms.
Furthermore, we also compared the PSTH of the spike trains generated by the OTP/BSG cascade with the PSTH of an Or59b OSN obtained from electrophysiology recordings in response to acetone concentration waveforms (Kim et al., 2015). As shown in Figure 4, the OTP/BSG closely matches the odorant response of the Or59b OSN.
Response of the Or59b OSN Model to White Noise Acetone Waveforms To further compare the response of the OTP/BSG cascade with the Or59b OSN response, we stimulated OTP/BSG cascades with white noise stimuli, and compared the PSTH of the model with the one from experimental recordings. The white noise stimulus was previous used in the experimental setting of (Kim et al., 2011b) for characterizing the response of Or59b OSNs to acetone. The output of each of the stages of the Or59b OSN model are shown in Figure 5.A. The odorant onset at around 1 second is picked up by the odorant concentration profile (see the the second row of Figure 5.A. In addition, the white noise waveform between 2 and 10 second is smoothed out. The smoothing effect is due to the peri-receptor process filter. The OTP model further emphasizes the gradient encoding (the third row of Figure 5.A), and predominantly defines the temporal response of the OSN model to white noise stimuli. The BSG output follows the OTP output, as the BSG is simply a sampling device. Lastly, we compare the model output and the PSTH from OSN recordings in (Kim et al., 2011b) (the fifth and the sixth rows of Figure 5.A). The Or59b OSN model output PSTH closely matches the PSTH obtained from recordings.
The peri-receptor process filter is critical in processing the white noise waveforms, but less critical in processing the static waveforms discussed in this work. The filter prevents the model from overemphasizing the gradient of the white noise waveforms. In absence of this filter, the response of the OTP/BSG cascade is severely limited in matching the response of Or59b OSN (Kim et al., 2011b) to acetone waveforms.
Response of the Or59b OSN Model to Staircase Acetone Waveforms Next, we stimulated OTP/BSG cascade with the staircase waveform that was previously used in experiments (Kim et al., 2011b), evaluated the PSTH from the resultant spike sequences, and compared the model PSTH to the one from experimental recordings.
As shown in the second row of Figure 5.B, the filter ℎ( ) in Equation 2 has negligible effect on the odorant concentration profile since the staircase is smooth unlike the white noise stimulus discussed above. The encoding at jump times is strongly sharpened by the OTP. Overall, the fruit fly OTP/BSG cascade indeed encodes both the concentration and concentration gradient. In particular, at each upward concentration jump, the PSTH of the OSN launches upon a local maximum and then drops down to a saturation point. In addition, at each downward concentration jump, the same PSTH drops down first to a local minimum and then bounces back.
In short, the OSN model closely reproduces the temporal response of Or59b OSNs for all tested stimuli. This suggests that the OPT/BSG cascade has the desired complexity to effectively model the fruit fly OSNs.

Evaluating Affinity, Binding and Dissociation Rates of Other (Odorant, Receptor) Pairs
We next interrogate the role of the binding and dissociation rates in the OTP/BSG cascade. For a given receptor type and two odorants with different binding rates and the same dissociation rate, responses of the OTP/BSG cascade are identical if the waveforms of two odorants only differ by a scaling factor that is the reciprocal of the ratio of the biding rates. This follows from Equation 3.

Response of Or59b OSNs to Two Different Odorants
We verify the prediction mentioned above by stimulating the Or59b OSN model with two odorant stimuli, acetone and 2-butanone, paired with four concentration waveforms, and compare the responses with the experimental recordings in (Kim et al., 2011a). As shown in the first row of Figure 6, the two odorants have identical normalized waveforms scaled by two different factors, 100 and 10. The affinity of acetone and 2-butanone were estimated to be 3.141 ⋅ 10 −4 and 3.164 ⋅ 10 −3 , respectively, and the dissociation rate of the two odorants were estimated to be 1.205 ⋅ 10 1 and 1.203 ⋅ 10 1 , respectively.
The two odorant stimuli elicit almost exactly the same response from the Or59b OSN recodings (see the second row of Figure 6) as well as from the output of the OTP/BSG cascade (see the third row of Figure 6). The difference in binding rate for acetone and2-butanone is perfectly counterbalanced by the scaling factors of the odorant waveforms in Figure 6. In addition, the output of the OTP/BSG cascade closely reproduces the PSTH of the Or59b OSN as shown in the forth row of Figure 6.

Evaluating the Odorant-Receptor Response of the OTP/BSG Cascade
We further investigated the role of the binding rate using three odorant-receptor pairs that were previously used in experimental settings (Kim et al., 2011a). In addition to the binding rate estimated for (acetone , Or59b), we applied Algorithm 1 to two additional odorant-receptor pairs using the original raw data presented in (Kim et al., 2011a): 1. for (methyl butyrate , Or59b) an affinity value 4.264 ⋅ 10 −4 and a dissociation value 3.788 ⋅ 10 0 were obtained from the steady-state spike rate at 87 spikes per second and the peak spike rate at 197 spikes per second in response to a constant stimulus with amplitude 20 ppm; 2. for (butyraldehyde , Or7a) an affinity value of 7.649 ⋅ 10 −3 and a dissociation value 8.509 ⋅ 10 0 were obtained from the steady-state spike rate at 43 spikes per second and the peak spike rate at 101 spikes per second in response to a constant stimulus with amplitude 173 ppm; We simulated the OSN model for each of the three odorant-receptor pairs with three types of stimuli, step, ramp, and parabola. The binding and dissociation rates for different odorant-receptor pairs above were separately set.
As shown in Figure 7, with only the change in the value of the binding and dissociation rates, the OSN model closely matches the OSN's response for all three tested odorant-receptor pairs. The results in Figure 7 suggests that a pair of binding and dissociation rates is capable of closely matching the temporal response of different temporal odorant concentration waveforms.

Estimating the Odorant-Receptor Affinity Matrix with DoOR Datasets
The DoOR database integrates OSN recordings obtained with different measurement techniques (Galizia et al., 2010;Münch and Galizia, 2016), including in situ spike counts (De Bruyne et al., 1999;Hallem et al., 2004;Smart et al., 2008) and calcium activity (Pelz et al., 2006), among others. Spike counts are directly available from OSN spike train recordings. Relating calcium activity to spike activity is, however, error prone. We consequently focus here on the odorant-OSN response datasets of the DoOR database that contain spike count information (Hallem and Carlson, 2006). These datasets currently contain spike counts of 24 OSN groups in response to 110 odorants with a constant amplitude of 100 pm. The spike count is color coded and depicted Figure 8.A. By employing Algorithm 1, we empirically estimated the affinity value for all 110 ⋅ 24 = 2, 640 odorant-receptors pairs. The estimated affinity corresponding to each entry of the spike rate matrix shown in Figure 8.A is depicted in Figure 8.B.
In summary, the binding and dissociation rate model together with the rest of the OTP/BSG cascade define a family of OSN models, and provide the scaffolding for studying the neural coding for odorant identity and odorant concentration in temporal domain at the OSN antennae population level.

Evaluating the Temporal Response of the AMP LPU
We investigated the temporal response of the OTP/BSG cascade to various odorant waveforms, including step, ramp, parabola, staircase, and white noise waveforms. In addition, we biologically validated the cascade with electrophysiological recordings of OSNs by demonstrating that the cascade is capable of reproducing the complex temporal responses of OSNs for multiple odorant receptor pairs.
Here, we study the temporal response of the AMP LPU, that consists of 50 OSN groups (Masse et al., 2009;Su et al., 2009). Each of 50 OSN groups consists of 25 OTP/BSG cascades (neurons) that express an unique receptor type. We tested the AMP LPU with the same staircase waveform as in Figure 5.B. For an assumed odorant, we assigned the same odorant-receptor affinity to OTPs in the same OSN group. The value of the affinity for each of the 25 OTP ranges between 2 ⋅ 10 −4 and 10 −2 with a step size of 2 ⋅ 10 −4 . The dissociation rate for all OTP models was set to 10 2 . For simplicity, we used the same set of parameters for all cascades across all OSN groups. The parameters of the OTP model are given in Table 2,  (Hallem and Carlson, 2006). Each column represents an odorant, and each row represents an OSN receptor type. (B) Each entry of the affinity matrix is estimated from each entry of the spike rate matrix using the inverse of the function empirically determined with Algorithm 1. Note the log-scale color map for the affinity values. and the parameters of the BSG model are listed in Appendix 2. From the spike sequences generated by the 25 cascades we evaluated the PSTH for each of the 50 OSN groups. Figure 9. Preview of the animation demonstrating the AMP LPU in response to a staircase concentration waveform. The animation was rendered by NeuroGFX (Ukani et al., 2019). Each of the OSN groups consists of 25 fruit fly OTP/BSG cascades. The PSTH for each of the OSN groups was evaluated from the spike sequences generated by 25 cascades. The affinity for each of the 50 OSN groups was assumed to be ranging between 2 ⋅ 10 −4 and 10 −2 with a step size 2 ⋅ 10 −4 . The dissociation rate for all OTP models was set to 10 2 . The rest of parameters of both OTP and BSG are given in Table 2 and Appendix 2, respectively. We visualize the 50 PSTHs and provide the preview and the link to the animation in Figure 9. The animation is rendered by NeuroGFX, a key component of FFBO (Ukani et al., 2019) (both NeuroGFX and FFBO are briefly presented in the Methods and Materials section). As shown, the top plot in the animation (and in Figure 9) shows the staircase odorant waveform, and the bottom plot shows the 3D view of the 50 PSTHs.
The resultant PSTHs exhibit distinct temporal responses across different OSN groups. Both the concentration and concentration gradient of odorants with (overall) high binding rate values is 2D encoded.
For receptors with extremely low (overall) binding rate values, the OTP/BSG cascade generates an output only after the concentration amplitude exceeds a certain value. For example, as shown in Figure 9, OSNs expressing receptors marked with orange color remain silent in the time interval between 2 to 6 seconds under a weak amplitude concentration stimulation. They start reacting to the odorant stimulus after 8 seconds as the amplitude increases from 80 ppm to 100 ppm. This closely matches the experimental recordings (Hallem and Carlson, 2006).
To further evaluate the AMP LPU, we used the affinity matrix estimated from the DoOR database (see Figure 8), and simulated 24 OSN groups in response to 110 different odorants. The dissociation rate for OSN groups was assumed to be 132, as it can not be estimated from the available records in the DoOR database. We applied the same staircase odorant waveform as above, and visualized the PSTH of OSN groups with an animation. In Figure 10, we provide the preview and the link to the animation. As shown, the top of the animation in Figure 10 shows the staircase odorant waveform as a function of time, and the bottom of the animation (also in Figure 10) shows the spike rate matrix for 24 OSN groups and 110 odorants at each time point. Each row of the matrix represents an OSN group, and each column of the matrix corresponds to an odorant. The animation demonstrates that the intensity of the spike rate matrix increases dramatically at each jump of the staircase waveform and drops down to a steady state value afterwards. This striking feature is due to the large value of the concentration gradient at transition times and clearly stands out in the animation. Figure 10. Preview of the animation of the spike rate matrix of 24 OSN groups in response to 110 odorants. Each of the OSN groups consists of 25 OTP/BSG cascades. The PSTH for each of the OSN groups is evaluated from the spike sequences generated by the 25 cascades. Each row of the matrix represents an OSN group, and each column of the matrix corresponds to an odorant. The affinity for each pairs of OSN groups and odorant is estimated using the DoOR database (see Figure 8). The dissociation of all OSN groups is assumed to be 132. (top) Staircase odorant waveform. (bottom) Dynamics of the spike rate matrix across time.

Discussion
Successful modeling of encoding of odorants by olfactory sensory neurons spread across the antennae and maxillary palps requires means of easily constructing and testing a range of hypotheses regarding the transduction of odorants into spike trains. The essential functionality of olfactory sensory neurons that we focussed on here is their concurrent encoding of both odorant identity and odorant concentration. To address these two functional aspects we presented an in-depth description of OSNs modeled as two stage samplers that quantitatively encode both the odorant identity and its concentration profile.
We devised a class of modular OSN models as a cascade consisting of an odorant transduction process and a biophysical spike generator. The OTP model consists of an active receptor model and a co-receptor channel model. The the BSG model employed here is based on the Connor-Stevens neuron model. After developing the OTP and BSG models, our focus was on the biological validation of the OTP/BSG cascade. To validate our modeling approach, we examined the response of the fruit fly OSN model to odorant waveforms that were previously used in experiments with different odorants and receptors, and compared the model responses with electrophysiology recordings. Our results show that the OTP/BSG cascade model proposed here indeed closely matches the complex temporal response of OSNs.

Limitation of the DoOR Database
The odorant-receptor affinity matrix together with the spike rate matrix provide the macroscopic I/O characterization of the OSN model at the population level. However, in the absence of additional information, such as the slope, width, or peak of the OSN response to the odorant onset, the dissociation rate can not be estimated with Algorithm 1 as such information is currently not available in the DoOR database. Thus, whereas previously the affinity and the dissociation rate are both estimated for multiple odorant-receptor pairs, only the affinity can be estimated for the odorant-receptor pairs recorded in the DoOR database. The dissociation rate together with the affinity are both required for reproducing the temporal response of OSNs. The latter alone can only characterize the steady state response. This illustrates some of the limitations of the DoOR datasets for characterizing the temporal response properties of OSNs, despite their richness for characterizing the steady state response of odorant-receptor pairs.

Combinatorial Coding of Odorant Identity is Concentration Dependent
The output of simulated OSN population demonstrates that the odorant identity is encoded in the set of odorant-activated OSN groups expressing the same receptor type. Different odorants evoke different sets of OSN groups, as shown in Figure 10. More importantly, the size of the set expands or reduces as the odorant concentration increases or decreases.

Complexity of the Input Space of Olfactory Neuropils
Our approach models on the molecular level the combinatorial complexity of the transformation taking place in Drosophila antennae OSNs. The transformation maps the odorant-receptor binding rate tensor modulated by the odorant concentration profile and the odorant-receptor dissociation rate tensor into OSN spike trains, respectively. The resulting concentrationdependent combinatorial code determines the complexity of the input space driving olfactory processing in the downstream neuropils, such as odorant recognition (Badel et al., 2016) and olfactory associative learning (Lin et al., 2014;Hige et al., 2015).

Future Work
The odorant receptor and the pheromone receptor share similar temporal variability in response to input stimuli (Gu et al., 2009), despite the differences in protein structure and chemical signaling between the two receptor families. Therefore, the fruit fly OTP model can be extended to model pheromone receptors. Pheromones can be modeled as odorants with an extremely high single receptor binding and dissociation rates.
The active receptor model can be readily extended to odorant mixtures. One interesting question is to study the odorant encoding of OSNs in the presence of a background odorant. Another potential direction is to investigate the "cocktail party" problem of odorant mixtures (Rokni et al., 2014). 1.570 ⋅ 10 1 cutoff frequency of the filter modeling the peri-receptor process 1 8.000 ⋅ 10 −1 slope of the transition region of the peri-receptor process filter 1.750 ⋅ 10 −1 scaling factor of the filtered odorant concentration gradient 2 8.877 ⋅ 10 1 rate of activation of the gating variable of the co-receptor channel 2 9.789 ⋅ 10 1 rate of deactivation of the gating variable of the co-receptor channel 3 2.100 ⋅ 10 0 rate of increase of the state variable of the calcium channel 3 1.200 ⋅ 10 0 rate of decrease of the state variable of the calcium channel 7.089 ⋅ 10 3 feedback strength from the calcium channel to the co-receptor channel 7.534 ⋅ 10 −2 value achieving the half-activation of the co-receptor channel 1 the Hill coefficient of the co-receptor channel 7.774 ⋅ 10 1 maximum current amplitude generated by the co-receptor channel

Methods and Materials I/O Characterization of the OTP/BSG Cascade
Is the OTP model given in Equation 7 capable of qualitatively reproducing transduction currents as those recorded in voltage clamp experiments? We empirically explore this question below. We first empirically tuned the parameters of the odorant transduction process model so as to generate similar transduction currents as recorded in the voltage-clamp setup published in (Cao et al., 2016). The binding rates and the dissociation rates of all OTP models were set to 1 and 132, respectively, and values of the other parameters are listed in Table 2.
We evaluated the model using step stimuli ( ), ramp stimuli ( ), and parabola stimuli ( ), chosen as, (1 − 5( − 2.3)), 2.3 ≤ < 2.5 0, ℎ , ( 1 1.9 ( − 0.5)) 2 , 0.5 ≤ < 2.4 (1 − 10( − 2.4)) 2 , 2.4 ≤ < 2.5 0, ℎ , where is a scalar ranging between 1 and 101 with a step size of 5. The response at the output of the peri-receptor process * ℎ, the odorant concentration profile [ ] , and the ratio of bound receptor [ 1 ] are shown in Figure 11. The slope of the rising phase of * ℎ after the onset of odorant is due to the effect of the filter ℎ( ). The odorant concentration profile [ ] encodes the gradient of the concentration for the step stimuli (see the chair-shaped response), but less so for the ramp and parabola stimuli. Lastly, the bound receptor [ 1 ] transforms the odorant concentration profile and maps it into a bounded range between 0 and 1.
As shown in Figure 11.A, the OTP model exhibits temporal response dynamics akin to the adaption phenomena reported in (Cao et al., 2016). Furthermore, we tested the OTP model with ramp and parabola stimuli of different concentration amplitude values, as shown in Figure 11. Compared with the stimulus response of the active receptor model, the response of the OTP model exhibits a complex temporal variability, that is sensitive to both the amplitude and the gradient of the odorant stimulus waveform. For example, as shown in Figure 11.B, the response of the OTP model to the ramp stimulus first increases linearly as the ramp stimulus increases, but then plateaus and remains constant as the gradient of the ramp stimulus is a constant. In addition, as shown in Figure 11.C, the response of the OTP model to the parabola stimulus roughly resembles a ramp function that closely matches with the gradient of the parabola stimulus.
The complex temporal response of the OTP model is due to the feedback received by the co-receptor channel from the calcium channel. Without the calcium channel feedback, the OTP model is reduced to a three-stage (peri-receptor processing, bound-receptor generator, and co-receptor channel) feedforward model. The feedback enables the OTP model to encode the odorant concentration profile components, i.e., both the filtered odorant concentration and concentration gradient. In addition, the nonlinearities embedded in the current generation of the co-receptor channel (see also Equation 6) acts as a normalization block, that facilitates the OTP model to map a stimulus with a wide range of amplitude values into a bounded transduction current.
The temporal response variability of the OTP/BSG cascade is similar to the transduction current generated by the OTP model. The similarity between the responses of the OTP model and the OTP/BSG cascade suggests that the temporal variability of the odorant concentration profile is primarily encoded in the OTP model. The BSG model is simply a sampling device mapping input current waveforms into spike trains.
Evaluating the Steady State Response of the OTP/BSG Cascade Figure 12. The transformation of the odorant concentration amplitude into steady-state spike rate by the OTP/BSG cascade for fixed values of the ligand-receptor affinity in response to 5-second-long constant stimuli. The parameters of the OTP model are given in Table 2, and the parameters of the BSG model are listed in Appendix 2. The amplitude of the constant odorant stimuli ranges between 10 −3 and 10 3 with a step size of 0.1 on the logarithmic scale. The spike rate is calculated in a window between 4 and 5 seconds.
We note that Equation 3 can be written as, . (11) where [ ] ∕[ ] is the odorant-receptor or ligand-receptor "affinity" (Kastritis and Bonvin, 2013). The active receptor model postulated in Equation 11 implies that in steady state the product between the odorant-receptor affinity and the odorant concentration profile is the main figure of merit for I/O characterization of the fruit fly OTP/BSG cascade. To study its mapping into spike rate, we simulated OTP/BSG cascades with constant stimuli, and evaluated the spike rate at steady state. The amplitude of step stimuli ranges between 10 −1 and 10 5 with a step size of 0.1 on the logarithmic scale. The affinity ranges between 10 −2 and 10 1 with a step size of 0.01 on the logarithmic scale. The parameters of the OTP model are given in Table 2, and the parameters of the BSG model are listed in Appendix 2. The step stimulus is 5 second long, and the OTP/BSG cascades reach steady state roughly after 3 seconds. We calculated the spike rate using a window between 4 to 5 seconds, and plotted the results in 2D in Figure 12. Note that the -axis in Figure 12 is on the logarithmic scale. As shown in Figure 12, for different values of the odorant-receptor affinity, the mapping of the concentration amplitude into spike rate shifts along the x-axis. A low affinity value requires a higher concentration amplitude value in order to elicit spikes above the spontaneous activity rate.
As shown in Figure 12, the transformation of the product between the odorant-receptor affinity and the concentration amplitude into spike rate resembles a sigmoidal function. The OTP/BSG cascade starts spiking only after the product exceeds a certain threshold value. For odorant-receptor pairs with a low affinity, the firing activity requires a larger minimal amplitude of concentration than for those with a higher affinity value. This, again, coincides with experimental findings that odorant-receptor pairs with lower affinity require higher odorant concentration values in order to elicit spiking activity (Hallem and Carlson, 2006).

Reproducing the 2D Encoding of the OSNs
To examine whether the fruit fly OTP/BSG cascade exhibits the 2D encoding property, we stimulated the cascade with the set of 110 triangular concentration waveforms that were previously used in experiments (Kim et al., 2011b) with Or59b and acetone. The triangular waveforms and their trajectories are plotted in Figure 13.A and Figure 13.C. We applied each of the triangular waveforms to 25 OTP/BSG cascades, and evaluated the PSTH using the spike train of all 25 cascades with a 20 ms bin size and 10 ms time shift between consecutive bins. The binding and dissociation rates of all OTP/BSG cascades was set to 3.141 ⋅ 10 4 and 1.205 ⋅ 10 1 , respectively. The parameters of the OTP model are listed in Table 2, and the parameters of the BSG model are listed in Appendix 2.
The responses of the OTP/BSG cascade are given in Figure 13. The PSTH of the OTP/BSG cascade in response to different waveforms is color-coded in both the 2D and 3D view, as shown in Figure 13.B and Figure 13.D, respectively. In addition, we applied the 2D ridge regression algorithm to identify a 2D encoding manifold that best fits the PSTHs. The manifold and its contour are depicted in Figure 13.F and Figure 13.E, respectively. Similarly to the case of the staircase waveform, the OTP/BSG cascade firing rate increases dramatically as the concentration increases.
As shown in Figure 13.F, a 2D encoding manifold in a concentration and concentration gradient space provides a quantitative description of the OTP/BSG cascade. Examining Figure 13.F, we note that the 2D encoding manifold is highly nonlinear and that the OTP/BSG cascade clearly encodes the odorant concentration and its rate of change. The OTP/BSG cascade responds very strongly to even the smallest positive values of the gradient and encodes only positive concentration gradients at low odorant concentrations. At high concentrations the OSN mostly encodes the odorant concentration. An algorithm to estimate the values of the odorant-receptor biding and dissociation rates may, therefore, • estimate the ligand-receptor affinity in steady state when the LHS of Equation 11 is zero for all values of the dissociation rate [ ] , and • estimate of the dissociation rate [ ] during a concentration jump assuming the value of the ligand-receptor affinity to be the one obtained in 1. above.
We describe the procedure above in more detail in Algorithm 1.
Algorithm 1 Estimation of the Affinity, Binding and Dissociation Rates 1: procedure (given step stimulus, steady state spike rate, and peak spike rate).

Simulation Setup in the Fruit Fly Brain Observatory Platform
The Fruit Fly Brain Observatory (FFBO) (Ukani et al., 2019) is an open-source platform for the emulation and biological validation of fruit fly brain models in health and disease. It provides users highly intuitive tools to execute neural circuit models on modern commodity hardware (Givon and Lazar, 2016). NeuroGFX is the key FFBO component supporting the implementation of AMP LPU. NeuroGFX conjoins the simultaneous operations of model exploration and model execution with a unified graphical web interface that renders interactive simulation results.
We evaluated the response to each of the stimuli by 50 neuron groups, each group consisting of the same OSN type. Each of 50 groups consisted of 25 OSNs, and in total there were 1, 250 OSNs. We then computed the PSTH for each of OSN groups using the resultant 25 spike sequences in each of the groups. The PSTH had a 20 ms bin size and was shifted by a 10 ms time interval between consecutive bins. The parameters of all OTP models are given in Table 2. The binding rate was separately set for each odorant-receptor pair. We used the same set of parameters for all 1250 cascades, but generated different sample paths for the Brownian motion term in Equation 12. The parameters of the BSG model are listed in Appendix 2. .
Lastly, the current generated by each of the channels is given by,