Interaction of cellular and network mechanisms for efficient - Loria

Document Sample
Interaction of cellular and network mechanisms for efficient - Loria Powered By Docstoc
					Interaction of cellular and network mechanisms for
efficient pheromone coding in moths
Hana Belmabrouka, Thomas Nowotnyb, Jean-Pierre Rosparsc, and Dominique Martineza,1
  Laboratoire Lorrain de Recherche en Informatique et ses Applications (LORIA) Unité Mixte de Recherche 7503, Centre National de la Recherche Scientifique,
54506 Vandoeuvre-lès-Nancy, France; bCentre for Computational Neuroscience and Robotics (CCNR), University of Sussex, Falmer, Brighton BN1 9QJ, United
Kingdom ; and cUnité Mixte de Recherche 1272 Physiologie de l’Insecte Signalisation et Communication (PISC), Institut National de la Recherche
Agronomique, 78026 Versailles, France

Edited by John G. Hildebrand, University of Arizona, Tucson, AZ, and approved October 27, 2011 (received for review July 29, 2011)

Sensory systems, both in the living and in machines, have to be                and evolve in time in a stimulus-specific manner (10–12). Using
optimized with respect to their environmental conditions. The                  time in such a manner as a coding dimension, however, poses
pheromone subsystem of the olfactory system of moths is a partic-              a serious problem for a flying insect. In natural plumes, odor
ularly well-defined example in which rapid variations of odor                   contacts with the antennae are brief and intermittent, with up to
content in turbulent plumes require fast, concentration-invariant              five contacts per 1 s and each contact lasting down to under 20 ms
neural representations. It is not clear how cellular and network               (13, 14). Consequently, responses in the olfactory system have to
mechanisms in the moth antennal lobe contribute to coding                      be fast. In Drosophila, for example, upwind surges generally begin
efficiency. Using computational modeling, we show that intrinsic                85 ms after ORN onset (15). During this time frame, the PNs
potassium currents (IA and ISK) in projection neurons may combine              cannot fire many more than 10 spikes (assuming mean rates of 120
with extrinsic inhibition from local interneurons to implement a dual          spikes/s and 0 ms latency). Recordings from moths in natural
latency code for both pheromone identity and intensity. The mean               conditions revealed that MGC-PNs produce 3–28 spikes after each
latency reflects stimulus intensity, whereas latency differences
                                                                               pheromone contact on the antennae, depending on stimulus in-
carry concentration-invariant information about stimulus identity.
                                                                               tensity (figures 3 and 4 in ref. 13). These observations put strong
                                                                               time constraints on the optimization of pheromone representa-
In accordance with physiological results, the projection neurons ex-
                                                                               tions for concentration-invariant recognition. Given the discon-
hibit a multiphasic response of inhibition–excitation–inhibition. To-
                                                                               tinuous nature of the odor cues, how does the MGC achieve stable
gether with synaptic inhibition, intrinsic currents IA and ISK account
                                                                               representations in less than 100 ms? How is the pheromone
for the first and second inhibitory phases and contribute to a rapid
                                                                               encoded with only a few spikes per PN?
encoding of pheromone information. The first inhibition plays the                  In this work, we address this fundamental question for the ol-
role of a reset to limit variability in the time to first spike. The second     factory system of the moth Manduca sexta, which is one of the best-
inhibition prevents responses of excessive duration to allow track-            characterized olfactory systems. In natural conditions, male
ing of intermittent stimuli.                                                   M. sexta moths are attracted by a blend of two main compounds of
                                                                               the sex pheromone in a proportion of 2:1 (16). The primary natural
olfaction pheromone plume     | odor coding | concentration invariance |       component is E,Z-10,12-hexadecadienal [bombykal (BAL)]. The
small conductance channel                                                      second natural component, E,E,Z-10,12,14-hexadecatrienal, is
                                                                               often replaced by a more stable chemical E,Z-11,13-pentadeca-
                                                                               dienal (C15) in laboratory experiments (17). For simplicity, we will
T   o attract conspecific males, female moths release volatile
    blends of a few chemical compounds, so-called sex phero-
mones. Most pheromone compounds are hydrocarbon molecules
                                                                               refer to both BAL and C15 as natural pheromone components.
                                                                               MGC-PNs are classified as specialists or generalists according to
emitted in the right proportions for biological relevance. Sympat-             their responses to single compounds (18, 19). Generalist PNs are
ric species often share the same chemicals but in different pro-               excited by BAL and C15, whereas specialist PNs are excited by one
portions [e.g., Heliothis zea and H. virescens (1, 2) or Yponomeuta            compound and inhibited by the other (19, 20). BAL or C15 spe-
cagnagellus, Y. irrorellus, and Y. plumbellus (3)]. To prevent cross-          cialist PNs fire synchronously, and their synchronization is sharp-
attraction, the pheromone blend is distinguishable not only by the             ened by the presence of both components in the pheromone
identity of the individual components but also by their precise                mixture (21). With the blend, they exhibit a multiphasic discharge
ratio. Behavioral experiments revealed that male moths are tuned               pattern, either inhibition–excitation–inhibition (I1/E/I2) or excita-
to the specific proportions emitted by their conspecific females.                tion–inhibition (E/I2) (18). The goal of our work was to highlight
Slight changes in component ratios, for example, can inhibit the               extrinsic and intrinsic neuronal mechanisms that would lead to the
anemotactic response of males of one species, while attracting                 experimentally observed distribution of response properties and to
males from another species (4–6). Pheromones are passively                     investigate how pheromone blends are encoded or represented by
transported in the air and mixed by atmospheric turbulence, and                such PN multiphasic discharge patterns. The origin of this multi-
therefore, during flight, male moths encounter pheromone fila-                   phasic patterning is unclear. I2, for example, was abolished by
ments with all components present in the same proportions as re-               certain GABAA antagonists [e.g., bicuculline (BIC)] and not by
leased (7) but over a broad range of concentrations (8, 9). In the             others [e.g., picrotoxin (PTX)] (22). Differences in BIC and PTX
race for mating, flying male moths have to solve a difficult pattern             sensitivity might be explained by the fact that a subunit of
recognition problem (i.e., recognize, in real time, pheromone                  the GABAA receptor is PTX-sensitive and BIC-insensitive (23).
compounds and their proportions independently of blend con-                    However, I2 was relatively unaffected by changes in chloride
centration). To do so, the information delivered by olfactory re-
ceptor neurons (ORNs) is integrated in the antennal lobe by a
neural network involving inhibitory local neurons (LNs) and ex-                Author contributions: H.B. and D.M. designed research; H.B. performed research; T.N. and
citatory projection neurons (PNs). All synaptic connections be-                J.-P.R. contributed new reagents/analytic tools; H.B. and D.M. analyzed data; and D.M.
                                                                               wrote the paper.
tween these neurons take place in a subset of enlarged, sexually
dimorphic glomeruli, the macroglomerular complex (MGC). The                    The authors declare no conflict of interest.
circuitry of the MGC is similar to the one of the generalist olfac-            This article is a PNAS Direct Submission.
tory subsystem. In the generalist subsystem, prolonged odor                    1
                                                                               To whom correspondence should be addressed. E-mail:
stimulations (up to several seconds) create local field potential               This article contains supporting information online at
oscillations and transient PN assemblies that synchronize to them              1073/pnas.1112367108/-/DCSupplemental.

19790–19795 | PNAS | December 6, 2011 | vol. 108 | no. 49                                                 
concentrations, which makes the role of GABAergic synapses               A                          C                        D
questionable (24). Interestingly, BIC has been shown to block a
calcium-dependent potassium channel, the small conductance
(SK) channel (25, 26), and Ca2+-dependent K+ currents as well as
transient A-type K+ currents have been identified in the PNs of
M. sexta (27, 28). Questions previously raised find specific for-
mulations + here. Do intrinsic K+ currents (A-type or Ca2+-de-
pendent K ) combine with network mechanisms (inhibition from
LNs) to shape the multiphasic response? How do inhibitory and
excitatory phases contribute to efficient coding of pheromone? Do
generalist and specialist PNs play different roles in this coding? In
this work, we address these questions by means of computational          B
modeling using a detailed MGC model of M. sexta moths.
MGC Model Reproduces Physiological Firing Patterns. Experimental
studies have revealed that the MGC of male M. sexta moths
contains two types of inhibitory LNs (IIa and IIb) (29). Both LN
types are multiglomerular but differ in their pattern of connec-
tivity and their response to the pheromone. LNs-IIa have very
dense branches in the MGC and exhibit brief and phasic excitatory
responses (figure 17 in ref. 29). On the contrary, LNs-IIb possess        Fig. 1. Response patterns in the intact MGC network model. (A) Typical LN
a less dense arborization, and their response to the pheromone           responses (86 LNs-IIa and 68 LNs-IIb). The LNs-IIa responded to pheromone-
blend is a long-lasting tonic discharge (figure 19 in ref. 29). MGC-      like stimuli with a short burst of spikes. On the opposite end, the LNs-IIb
                                                                         exhibited prolonged tonic responses. (B) Typical PN response patterns
PNs are classified as specialist or generalist according to their
                                                                         extracted from the MGC model (41 PNs). From their response to single
responses to single compounds (18, 19). Generalist PNs are ex-
                                                                         compounds, the PNs were classified as specialist or generalist. Generalist PNs
cited by both compounds C15 and BAL (n = 13; 32% of the PNs in           (30% of the PNs; dark area in pie chart) were excited by both compounds,
table 1 in ref. 18). Specialist PNs are inhibited by one compound        C15 and BAL, whereas specialist PNs were inhibited by one compound and
and excited by the other. They exhibit a complex discharge pattern       excited by the other. When we looked at the network connectivity, we
with the blend of either excitation–inhibition (biphasic pattern E/I2    noticed that generalist and specialist PNs were multiglomerular and unig-
in 39% of the PNs) or inhibition–excitation–inhibition (triphasic        lomerular, respectively. Specialist PNs responded to the blend with a com-
pattern I1/E/I2 in 29% of the PNs).                                      plex temporal pattern, either biphasic with excitation–inhibition, denoted E/
   As described in Methods, we developed a detailed model of the         I2 (37% of the PNs; gray area in pie chart), or triphasic with inhibition–ex-
M. sexta MGC that fits the numbers used in refs. 18 and 29 (i.e., 86      citation–inhibition, denoted I1/E/I2 (33% of the PN; white area in pie chart).
LNs-IIa, 68 LNs-IIb, and 41 PNs). A schema of the model is shown         The proportions of response patterns, 30%–37%–33% for the generalists
in Fig. S1. We investigated whether this model reproduces the            and biphasic and triphasic specialists, were obtained on 30 MGC networks
different types of responses in the correct proportions. In line with    randomly generated with p(LN-IIa → PN) = 0.5 and p(LN-IIb → PN) = 0.2 and
the experimental data (figures 17 and 19 in ref. 29), LNs-IIa and         simulated during 1 s of biological time. The real proportions (areas delimited
LNs-IIb in our MGC model responded to blend-like stimuli, with           by red lines in the pie chart) observed in experimental data are 32%–39%–
short- and long-lasting firing patterns, respectively (Fig. 1A).          29% (29). The Kullback–Leibler distance (KLD) between the real proportions
                                                                         and those proportions obtained with the model was KLD = 1.33. (C) Capa-
Generalist and specialist responses of the PNs were also repro-

                                                                         bility for the MGC model to reproduce the physiological responses for dif-
duced (Fig. 1B). Their proportions, however, depended on the
                                                                         ferent probabilities of connection p(LN-IIa → PN) and p(LN-IIb → PN). The
probabilities of connection p(LN-IIa → PN) from LNs-IIa to PNs           KLD is zero when the model reproduces the physiological responses in the
and p(LN-IIb → PN) from LNs-IIb to PNs. To quantify the plau-            right proportions. The proportions of response patterns were obtained, for
sibility of the model, we computed the Kullback–Leibler distance         each pair of probabilities p(LN-IIa → PN) and p(LN-IIb → PN) on 30 random
(KLD) between the proportions of PN responses obtained with the          MGC networks (simulations of 1 s). (D) The MGC model is optimal for p(LN-
model and the real proportions observed in experimental data             IIa → PN) = 0.5 and p(LN-IIb → PN) = 0.2. At this point, KLD = 1.33, and the PN
(gray-coded areas vs. areas delimited by red lines in the pie chart in   responses with their proportions are illustrated in B.
Fig. 1B). The KLD equals zero when the model reproduces exactly
the biological responses in the right proportions. For each pair of
probabilities p(LN-IIa → PN) and p(LN-IIb → PN), 30 MGC                  LNs-IIb, we selectively removed each type of LNs from the MGC
networks were randomly generated and simulated during 1 s bi-            network and analyzed the PN responses. Without LN-IIa, the tri-
ological time. The PNs were classified automatically with an al-          phasic patterns were lost (Fig. S2A). The biphasic patterns rep-
gorithm identifying their responses to the blend and single              resented 68% of the PNs, which corresponds to the proportions of
compounds. The KLD (mean and SD) was estimated over the                  biphasic plus triphasic patterns in the intact network. These results
different runs (Methods and SI Text, Model). As expected, the KLD        indicate that removing LNs-IIa from the network changed the
depended on LN to PN connectivity (Fig. 1C). The existence of            patterns from triphasic (I1/E/I2) to biphasic (E/I2). LNs-IIa are,
a minimum at p(LN-IIa → PN) = 0.5 and p(LN-IIb → PN) = 0.2               therefore, responsible for I1 in the triphasic patterns. As for I2, the
suggests that the PNs should be more densely connected to the            analysis did not reveal any particular role for LNs-IIa. The same
LNs-IIa than the LNs-IIb to reproduce physiological firing patterns       procedure was applied to study the role of the LNs-IIb. Removing
(Fig. 1D). This result is in line with experimental data revealing       the LNs-IIb did not disturb the proportions of the PNs responses
a dense and a sparse glomerular innervation for LNs-IIa and LNs-         (same pie chart as for the intact network in Fig. 1B). However,
IIb, respectively (29). Although the minimal KLD is 1.33 and not 0,      although I2 did not disappear from the multiphasic responses, its
the PN responses and their proportions obtained in simulation are        duration was significantly reduced (106 ± 2.68 ms without LNs-IIb
in good agreement with those obtained experimentally (Fig. 1B).          vs. 112 ± 5 ms for the intact network, P < 0.05, Kruskal–Wallis
                                                                         test) (Fig. S2B). We, therefore, concluded that the LNs-IIb are
Roles of Extrinsic Properties: LNs-IIa and LNs-IIb. From the previous    involved in but not solely responsible for I2.
section, it becomes clear that inhibition from LNs plays a role in          Another putative role for the LNs is to synchronize the PNs. In
shaping the PN response. However, LNs-IIa and LNs-IIb have               the MGC of male M. sexta, BAL or C15 specialist PNs fire in syn-
different firing responses and connectivity patterns (Fig. 1 A and        chrony when activated by pheromone stimulation to the antenna
D). To separate the potentially different influences of LNs-IIa and       (21). Interestingly, greater synchrony is obtained with the blend

Belmabrouk et al.                                                                            PNAS | December 6, 2011 | vol. 108 | no. 49 | 19791
than single compounds. To assess whether our model can re-                                                           could explain the significant delay to the first spike in PNs com-
produce these data, we examined the precise timings of spikes and                                                    pared with LNs? When the A current was blocked, PNs fired
computed the level of synchrony as the percentage of synchronized                                                    earlier (latency was 54 ± 18 ms for intact PNs and 26 ± 22 ms with
spikes within 150 ms after stimulus onset (Methods). The data                                                        IA blocked, P < 0.05, Kruskal–Wallis test) (Fig. 2A). Thus, the
reported in Fig. S2C are means and SDs estimated over 30 trials. In                                                  firing delay in PNs resulted from the action of the A current. Fast
the intact MGC model, BAL or C15 PNs are more synchronized                                                           activation of the transient K+ current prevented the PNs from
with the blend (69 ± 27%) than with single compounds (59 ± 26%).                                                     firing immediately in response to a depolarization. When the A
Simulation data are in agreement with experiments (figure 4c in                                                       current was blocked, the PNs were also less synchronized (syn-
ref. 21). Is greater PN synchrony with the blend attributable to                                                     chrony level = 69 ± 27% for intact PNs and 66 ± 23% with IA
more vigorous inhibition from more strongly activated LNs? To                                                        blocked, P < 0.05, Kruskal–Wallis test) (Fig. 2B). The decrease in
answer this question, we selectively removed each type of LNs from                                                   synchrony when IA is blocked is, therefore, attributable to earlier
the MGC network. The consequence was a decrease in synchrony                                                         PN firings, which makes the reset by I1 less effective.
from 69 ± 27% for the intact MGC model to 51 ± 23% for the                                                              Previous simulations highlighted the causes and effects of I1. As
model without LN-IIa and 56 ± 28% for the model without LN-IIb                                                       for I2, experimental data revealed that it is abolished by BIC and
(Fig. S2C) (P < 0.05, Kruskal–Wallis test followed by Mann–                                                          not PTX (22) and that a calcium-dependent potassium channel,
Whitney test). Although both LN types are involved in PN syn-                                                        the SK channel, is BIC-sensitive (25, 26). Given that I2 remained
chrony, removing the LNs-IIa had significantly more impact. The                                                       relatively unchanged in the MGC model without LNs, as seen
LNs-IIa contribute more to PN synchrony, because the PNs are                                                         above, could it be produced by the SK current? Blocking ISK did
more densely connected to the LNs-IIa than to the LNs-IIb [p(LN-                                                     not disturb the proportions of the PN responses (same pie chart as
IIa → PN) = 0.5 vs. p(LN-IIb → PN) = 0.2]. To understand how                                                         for the intact network in Fig. 1B), but I2 disappeared from the
LNs-IIa synchronize the PNs, we analyzed the PN activities in the                                                    multiphasic responses (Fig. 2C). Triphasic (I1/E/I2) patterns were
intact MGC network and the same network with the LNs-IIa                                                             changed to biphasic (I1/E), and biphasic (E/I2) patterns were
blocked. In the intact network, we observed that inhibition I1 arose                                                 changed to monophasic ones (E). In all cases, the PN response
before the E response (representative voltage traces of two PNs are                                                  without the SK current consisted in a long tonic excitation (E
represented in Fig. S2D). During the I1 phase, the PNs have                                                          phases in Fig. 2C), which disrupted the ability to track pheromone-
a tendency to forget their initial conditions by relaxing their activity                                             like pulses (Fig. 2D). Given that the SK channel is sensitive to BIC,
to a hyperpolarized quasi-steady state. As a consequence, the                                                        these results are in agreement with experimental data, showing
neurons fire synchronously when inhibition is removed. On the                                                         that BIC prolongs PN excitation and compromises coding of in-
contrary, in the network where the LNs-IIa were blocked, the I1                                                      termittent pheromone pulses in M. sexta (30).
phase disappeared or was not strong enough to play the role of
a reset and eliminate the influence of initial conditions (Fig. S2D).                                                 Efficient Coding of Pheromone: Dual Latency Coding Scheme. It is
                                                                                                                     known that behavioral responses of male moths to pheromone
Roles of Intrinsic Properties: IA and ISK. The above observations in-                                                signals are fast. However, the nature of the neural code that
dicate that enhancement of PN synchrony is attributable to a loss                                                    underlies rapid pheromone processing has been elusive. To ex-
of initial conditions and can be interpreted in terms of a reset by                                                  amine candidate neural codes, we simulated the MGC model with
inhibition. A prerequisite to have synchronization in the first place                                                 inputs mimicking a mixture of BAL and C15 in different pro-
is that inhibition from LNs arises before the PN response (I1 be-                                                    portions and various concentrations. The binary mixture model
fore E), which implies that LNs need to fire before PNs. In                                                           was given by K(P × BAL + (1 − P) × C15), with K being the
agreement with this hypothesis, first spike latencies were signifi-                                                    concentration and P being the proportion (SI Text, Model). In our
cantly shorter for LNs than PNs (19 ± 3 ms for LN-IIa and 3 ± 1.5                                                    model, P = 2/3 corresponds to the correct natural pheromone
ms for LN-IIb vs. 54 ± 18 ms for PNs, P < 0.05, Kruskal–Wallis test                                                  ratio (proportion of BAL two times as big as C15). Spike raster
followed by Mann–Whitney test) (Fig. 2A). What intrinsic current                                                     activities revealed three distinct populations of synchronized

A                   80                                               B                      100
                                                                         Specialist PN−PN
                                                                          Synchrony (%)

                                                                                            80      a
     Latency (ms)

                    60                           c

                    40                                                                      40
                         a                                                                  20
                    20                                                                                                           Fig. 2. Effects of removing IA or ISK. (A) First spike latencies
                                   b                                                         0                                   were significantly shorter for LNs than PNs (19 ± 3 ms for LN-IIa
                    0                                                                              Intact   I   blocked
                         LN−II LN−IIb          PN         PN
                                                                                                            A                    and 3 ± 1.5 ms for LN-IIb vs. 54 ± 18 ms for PNs, P < 0.05,
                               a            (I Active)(IA blocked)
                                             A                                                                                   Kruskal–Wallis test followed by Mann–Whitney test). Without
C                                                                                                                                IA, the PNs fired earlier (latency = 26 ± 22ms, P < 0.05, Kruskal–
                                                                                                                                 Wallis test). (B) Without IA, PNs were less synchronized. The
                                                                                                                                 percentage of synchronous spikes found in all pairs of specific
                                                                                                                                 PNs was 69 ± 27% for intact PNs vs. 66 ± 23% without IA (P <
                                                                                                                                 0.05, Kruskal–Wallis test). (C) Without ISK, I2 disappeared from
                                                                                                                                 the multiphasic responses. The original biphasic and triphasic
                                                                                                                                 patterns in B were changed, respectively, to monophasic (E)
                                                                                                                                 and biphasic (I1/E). The same proportions of 30%–37%–33% as
                                                                                                                                 in B were, thus, obtained for the generalists and monophasic
                                                                                                                                 and biphasic specialists, respectively. (D) Without ISK, PNs lost
D                    1                           SK blocked                                                                      pulse-tracking capability. The autocorrelation function was
              0.8                                Intact network
                                                                                                                                 computed to quantify the capability of PNs to track phero-

              0.6                                                                                                                mone-like pulses delivered at 1.67 Hz. For intact PNs, the au-
                                                                                                                                 tocorrelation function (black curve) presented periodic peaks
                                                                                                                                 separated by 600-ms intervals and corresponding to the stim-
              0.2                                                                                                                ulus interpulse interval. Without ISK, the autocorrelation
                     0                                                                                                           function was not periodic anymore (red curve). The alteration
                          −2000         0        2000
                             time lag (in msec)                                                                                  of the ability to track pheromone-like pulses is illustrated by
                                                                                                                                 comparing the response of a PN with or without ISK (Right).

19792 |                                                                                                                            Belmabrouk et al.
neurons: generalist PNs, specialist C15-PNs, and specialist BAL-                                                                                        concentration (red curve in Fig. 3A). Despite relative concen-
PNs (Fig. S3A). Generalist PNs were the most activated neurons,                                                                                         tration invariance, information on concentration was not lost.
because they received inputs from both glomeruli. As a conse-                                                                                           The mean absolute latency of the population of generalist PNs
quence, they fired first, whereas specialist PNs fired later. With                                                                                         was found to be sensitive to changes in concentration K (black
a mixture proportion of P = 2/3, the BAL-PNs were followed by the                                                                                       curve in Fig. 3A), while being invariant to changes in proportion
C15-PNs, because they received two times as much excitatory drive                                                                                       (black curve in Fig. 3B). The MGC could, therefore, make use of
as the C15-PNs. The rank order of the three PN populations was                                                                                          latencies to implement a dual coding scheme, quantitative and
consistent across repeated trials, despite random initial conditions                                                                                    qualitative (a theoretical explanation for this dual invariance
thanks to the resetting of the PNs by I1 (see above). In the context of                                                                                 coding scheme is provided in SI Text and Figs. S4–S6). However,
fast processing, we asked whether first spike latency of the different                                                                                   latency coding makes sense only if one considers the first spike
PN populations could be used as a code. We, therefore, evaluated                                                                                        with respect to a reference signal [e.g., the first spike after
the information contained in a latency code compared with a rate                                                                                        stimulus onset (absolute latency as above) or after a particular
code. The datasets for rate or latency contained three attributes                                                                                       event in the brain dynamics (relative latency)]. Because the
(mean absolute latencies or mean firing rates of generalist, specialist                                                                                  stimulus onset is generally unknown to the brain, we looked at
C15, and specialist BAL populations), and the goal was to dis-
                                                                                                                                                        a putative internal reference signal from which to extract relative
criminate the correct proportion of P = 2/3 against incorrect mix-
                                                                                                                                                        latencies. We found that the generalist PNs fired consistently
ture ratios (Methods). Both rate and latency datasets were almost
linearly separable, which was revealed by plotting the average dis-                                                                                     with a delay of ∼40–50 ms after stimulus onset. Compared with
criminant plane obtained by cross-validation (Fig. S3 B and C).                                                                                         the specialist PNs, this delay was relatively constant across a wide
Classification performance was estimated by averaging the classifi-                                                                                       range of proportions and concentrations (Figs. S5 and S6),
cation errors over the validation sets. It was higher with latency                                                                                      suggesting that the firing of the generalist PN population could
attributes than rate attributes, with an average 99.7% classification                                                                                    provide an internal reference for temporal marking. We, there-
success with latency compared with 96.4% with rate (Fig. S3D). This                                                                                     fore, computed relative latencies of C15 and BAL PNs with
result showed that a code based on latency would be efficient for                                                                                        respect to the generalist PN population. We found that the
rapid pheromone recognition, but it does not give any clue as to how                                                                                    normalized difference in relative latencies encoded for the
concentration-invariant recognition could be achieved.                                                                                                  stimulus identity, being concentration-invariant (red curve in Fig.
   Interestingly, the normalized latency difference between C15                                                                                         3C) and proportion-sensitive (red curve in Fig. 3D), whereas the
and BAL populations, computed as (LBAL − LC15)/(LBAL +                                                                                                  mean relative latency of specialist PNs encoded for the stimulus
LC15), was found to be sensitive to the proportion P (red curve                                                                                         intensity, being concentration-sensitive (black curve in Fig. 3C)
in Fig. 3B), behaving as 2P − 1, but changed very little with                                                                                           and proportion-invariant (black curve in Fig. 3D).

              A                                               Blend = K*[P*BAL+(1−P)*C15]                                                               B                                                                  Blend = K*[P*BAL + (1−P)*C15]
                                                  60                    (P = 2/3)                   1                                                                                                      60                                                               0.5
                Mean generalist PN latency (ms)

                                                                                                                                                                                                                                       (K = 2)
                                                                                                                                                                         Mean generalist PN latency (ms)

                                                                                                                                                                                                                                                                                          Normalized difference latency
                                                                                                                        Normalized difference latency

                                                                                                                                                                                                                              Absolute Latencies
                                                               Absolute Latencies

                                                  45                                                0                                                                                                      45                                                               0

                                                  30                                                −1                                                                                                     30                                                               −0.5
                                                        K=1         K=2         K=3           K=4                                                                                                                P = 0.3     P = 0.4    P = 0.5    P = 0.6    P = 0.7

              C                                               Blend = K*[P*BAL + (1−P)*C15]                                                             D                                                                  Blend = K* [P*BAL+(1−P)*C15]                 1
                                                                         (P = 2/3)                                                                                                                                                     (K = 2)
                                                                                                                                                                                                                                                                                Normalized difference latency

                                                  150                                               1
                                                                                                                                                        Mean relative latency (ms)
                                                                                                         Normalized difference latency
              Mean relative latency (ms)

                                                                   Relative Latencies                                                                                                                                      Relative Latencies

                                                                                                                                                                                         40                                                                             0
                                                  90                                                0

                                                                                                                                                                                                           0                                                            −1
                                                  30                                                −1                                                                                                          P = 0.3     P = 0.4    P = 0.5    P = 0.6    P = 0.7
                                                        K=1         K=2        K=3          K=4

Fig. 3. Response latencies of the MGC model to blend-like stimuli in different proportions P and concentrations K. The binary mixture model is K(P × BAL +
(1 − P) × C15). (A) Dependence of concentration K on latency. On the one hand, the latency of the generalist PNs multiglom was sensitive to changes in
concentration K (black curve). On the other hand, the normalized latency difference (LC15 − LBAL)/(LBAL + LC15) between C15 and BAL PNs changed very
little with concentration (red curve). (B) Dependence of proportion P on latency. The normalized latency difference behaved as 2P − 1 (red curve), whereas the
latency of the generalist PNs was invariant to changes in proportion (black curve). (C and D) Same as in A and B but for the relative latencies of C15 and BAL
PNs computed with respect to the generalist PN population.

Belmabrouk et al.                                                                                                                                                                                                             PNAS | December 6, 2011 | vol. 108 | no. 49 | 19793
Discussion                                                              forward, lacking excitatory connections from PNs to LNs. The
When flying in a turbulent plume, male moths encounter inter-            LNs are, thus, driven by ORN inputs and adapt rapidly. Con-
mittent pheromone patches. Their short duration puts a strong           sequently, the PNs are reset by a single strong inhibitory phase I1
temporal pressure on the olfactory system (13, 14). In the context of   and fire synchronously when I1 ends. This type of nonoscillatory
fast pheromone processing, we suggested that latency to the first        synchronization persisted with added excitatory connections
spike could be used as information carrier. Evidence for latency        from PNs to LNs, provided that they are weak enough (strong
coding has been found in diverse neuronal structures such as the        PN-LN excitatory connections elicited fast ∼40-Hz oscillations in
visual system (31–33), the tactile system (34), and the olfactory       our model) (SI Text and Fig. S7).
system (35–37). Latency coding often implies sensing the external          A prerequisite for the inhibitory reset I1 is that LNs fire before
world as a discrete sequence of events. In the visual system, sudden    PNs. In our MGC model, spike latencies were significantly shorter
movements of the eyes, called saccades, sample the visual input         for LNs than PNs. PN firing was delayed by about 50 ms with re-
into discrete snapshots. After each saccade, a new image is encoded     spect to the fastest LNs, and this delay was attributable to the A
into first spike latencies to be further processed by the brain.         current in PNs (Fig. 2A). In line with our results, time lags of about
Similarly to the saccades in vision, sniffing in rodents samples the     60 ms have been reported between LN and PN firing latencies in
olfactory scene into discrete odor puffs. In each of these cases, the   the bee antennal lobe (49). Transient A-type K+ currents have
stimulus onset is actively controlled by the brain. In the case of      been identified in the PNs of M. sexta moths (27, 28). However,
moths, however, the pheromone plume is broken up into discrete          characterizing their potential involvement in modulating the time
filaments by the turbulent medium and not by an internal sampling        to first spike needs additional experimental studies.
mechanism. Under this condition, what could be the internal ref-           The various response patterns observed experimentally in M.
erence signal from which to extract latencies? Our MGC model            sexta projection neurons were reproduced in our MGC model by a
suggests that generalist PNs serve as a temporal marker to compute      combination of intrinsic K+ currents (A-type or Ca2+-dependent
relative latencies of specialist PNs. Interestingly, normalized la-     K+) and network mechanisms (inhibition from LNs) rather than
tency differences of specialist PNs carried concentration-invariant     by varying the connectivity as in ref. 50. These intrinsic currents
information about pheromone blend identity. Subtraction and             allowed the PNs to encode fast pulses of pheromone, despite long-
normalization to achieve invariance were previously reported. In        lasting ORN responses. The ability to track intermittent stimuli
the salamander visual system, differences in response latency of        was largely because of an intrinsic Ca2+-dependent K+ channel
retina cells were found to be invariant to stimulus contrast (33). In   that produced I2 and thereby, prevented PN responses of excessive
the rat olfactory bulb, normalization of glomerular responses has       duration. Ca2+-dependent K+ currents have been found in the
been proposed as a means to achieve invariance in concentration         PNs of M. sexta moths, but their type was not identified (28). A
(38), and we have suggested a mechanism of competition by lateral       Ca2+-dependent K+ current has been incorporated in another
inhibition in the moth MGC in a separate model (39). From this          modeling study as a mechanism of spike frequency adaptation to
coding perspective, the novelty of the work presented here lies in      reproduce the slow patterning observed in the generalist system
the use of response latencies for signaling the correct component       (51). For M. sexta moths, our model suggests that the Ca2+-de-
ratio rather than the presence of network oscillations (40) or spe-     pendent K current involved in I2 could be of the SK type. This
cific rate patterns (41). A latency code allows much faster responses    prediction is supported by experimental data revealing that I2 is
and therefore, seems best suited to address the tight time con-         abolished by BIC and not PTX (22), a difference explained by the
straints of pheromone blend detection in a turbulent plume. An-         fact that the SK channel is BIC-sensitive (25, 26). To our knowl-
other important point in our model is that, despite relative            edge, SK channels in insects have been shown only in a visual
concentration invariance, information on concentration was not          interneuron of the locust (52). A direct confirmation of the exis-
lost as the mean latency of specialist PNs encoded for stimulus         tence of SK channels in M. sexta moths is, therefore, needed.
intensity. The MGC of M. sexta moths could, therefore, make use of      However, this task is not easy to accomplish, because the most
latencies to implement a dual coding scheme and resolve the ap-         common antagonist, apamin, has no effect in insects (52).
parent contradiction between quantitative coding (sensitive to
changes in concentration and invariant to changes in proportion)        Methods
and qualitative coding (sensitive to changes in proportion and in-      Network Simulations. The simulated MGC model consisted of 86 LNs-IIa, 68
variant to changes in concentration). As an extension to this work,     LNs-IIb, and 41 PNs, and these numbers were grounded experimentally (18,
the proposed latency coding scheme provided a direct input for          29). LNs and PNs were modeled as one compartment conductance-based
designing bio-inspired data analysis methods for artificial olfaction    models involving ionic currents and leak (SI Text, Model). In M. sexta, the LNs
in electronic noses (42).                                               are multiglomerular, but most of the PNs are uniglomerular (53). We,
   In latency coding, one considers the neurons as analog to delay      therefore, modeled the LNs as multiglomerular. Uni- and multiglomerular
converters: the most strongly activated neurons tend to fire first,       PNs were randomly chosen with probabilities of 0.7 and 0.3, respectively.
whereas more weakly activated cells fire later or not at all.            Uniglomerular PNs were connected with equal probability to either one or
Random initial conditions may, however, jeopardize this ideal           the other glomerulus (toroid-BAL or cumulus-C15). The glomerular input
situation by introducing variability in the time to first spike. The     model is detailed in SI Text, Model. Inhibitory LNs, IIa and IIb, were randomly
simulations presented in this article showed that this condition is     connected to the PNs with probabilities of p(LN-IIa → PN) and p(LN-IIb → PN).
not the case. Inhibition I1 from LNs before the PN responses            The effect of these probabilities on the firing response patterns was studied
                                                                        above (Results). Physiological firing patterns were best reproduced with the
played the role of a reset by transiently driving the PNs to a
                                                                        intact MGC network and p(LN-IIa → PN) = 0.5 and p(LN-IIb → PN) = 0.2. These
common hyperpolarized state. As a consequence, PN latencies
                                                                        probabilities were also used in the simulations with nonintact networks. The
after stimulus onset depended more on the strength of the input         GABAA synaptic model from LNs to PNs was adapted from ref. 54 and is
than initial conditions. Generalist PNs fired first followed by           detailed in SI Text, Model. Network simulations were performed with SIR-
BAL and C15 PNs (Fig. S3A). The rank order of the three                 ENE, a C-based neural simulator developed by our team and available at
groups reflected the intensity of their inputs. Within each group, The ordinary differential equations are numer-
PN synchrony arose as a consequence of inhibitory reset and             ically integrated with a classical fourth-order Runge–Kutta method with
shared ORN excitation. This type of synchronization is non-             a time step of 0.01 ms.
oscillatory and was observed experimentally in the MGC of M.
sexta moths (21, 43). It differs from the 30- to 40-Hz oscillatory      Data Analyses. We identified the PNs’ types from the response to single
synchronization found in the generalist olfactory system of the         compounds. PNs were considered as specialist or generalist according to
same moth species (44). With nonpheromonal odors, the LNs               whether they were excited by one compound only or by the two compounds
are mainly driven by the PNs, and therefore, PN-LN interplay            presented individually. Excitation E was detected if at least one spike was
generates rhythmic inhibition and creates field potential oscil-         found within 100 ms after stimulus onset. Inhibition I1 was detected
lations (45–48). In contrast, our MGC model is essentially feed-        whenever a hyperpolarization of the membrane potential occurred right

19794 |                                                                                 Belmabrouk et al.
after stimulus onset. To determine hyperpolarization, the membrane po-                        0.6, 0.7 and K = 1). The dataset for rate or latency contained 180 samples of
tential was filtered (fourth-order low-pass Butterworth, 50-Hz cutting fre-                    three inputs (mean absolute latencies or mean firing rates of generalist, spe-
quency) and averaged before and after stimulus onset (averaging window                        cialist C15, and specialist BAL populations). The training set consisted of 90
of 200 ms before and 80 ms after stimulus onset). Inhibition I2 was de-                       samples taken randomly from the entire dataset, and the test set consisted of
termined from visual inspection. To quantify the plausibility of the model,                   the remaining 90 samples. Training and test sets were resampled 10 times to
the proportion of response patterns (combination of I1, E, and I2 ) was                       estimate the cross-validation error. A linear support vector machine (SVM)
compared with the real proportion observed in experimental data by com-
                                                                                              with parameter one was used as classifier (55). All data analyses and statistical
puting the KLD (SI Text, Model). Synchrony between two PNs was quantified
                                                                                              tests were performed in Matlab ( except for the classi-
as the percentage of synchronous spikes (maximum time lag ± 2.5ms) rela-
                                                                                              fication performance, which was assessed with a custom-made SVM program.
tive to the total number of spikes as in ref. 21. The synchrony level was
expressed as the mean synchrony computed for all pairs of uniglomerular
PNs. Classification performance was estimated for the two coding schemes:                      ACKNOWLEDGMENTS. We are grateful to Dr. Philippe Lucas for helpful
                                                                                              discussions. This work was funded by the French and British national agencies
rate and latency. Absolute latency and mean firing rate were computed as,
                                                                                              [Agence Nationale de la Recherche (ANR), Biotechnology and Biological Sciences
respectively, the time to first spike and the number of spikes within a 200-ms                 Research Council (BBSRC)] in a coordinated manner for Collaborative Research
window after stimulus onset. The classification task was to discriminate the                   in Systems Biology Grant BSYS-006-02-PHEROSYS and European Research
correct proportion of P = 2/3 at concentrations K (80 samples for P = 2/3 and                 Project “Biologically inspired computation for chemical sensing” Grant FP7-
K = 1, 2, 3, 4) against incorrect mixture ratios (100 samples for P = 0.3, 0.4, 0.5,          216916-NEUROCHEM.

 1. Pope MM, Gaston LK, Baker TC (1982) Composition, quantification and periodicity of         27. Mercer AR, Hayashi JH, Hildebrand JG (1995) Modulatory effects of 5-hydroxytryp-
    sex pheromone gland volatiles from individual Heliothis virescens females. J Chem             tamine on voltage-activated currents in cultured antennal lobe neurones of the
    Ecol 8:1043–1055.                                                                             sphinx moth Manduca sexta. J Exp Biol 198:613–627.
 2. Pope MM, Gaston LK, Baker TC (1984) Composition, quantification, and periodicity of sex    28. Mercer AR, Hildebrand JG (2002) Developmental changes in the density of ionic
    pheromone volatiles from individual Heliothis zea females. J Insect Physiol 30:943–945.       currents in antennal-lobe neurons of the sphinx moth, Manduca sexta. J Neurophysiol
 3. Löfstedt C, Herrebout WM, Menken SBJ (1991) Sex pheromones and their potential                87:2664–2675.
    role in the evolution of reproductive isolation in small ermine moths (Yponomeuti-        29. Matsumoto SG, Hildebrand JG (1981) Olfactory mechanisms in the moth Manduca
    dae). Chemoecology 2:20–28.                                                                   sexta: Response characteristics and morphology of central neurons in the antennal
 4. Glover TJ, Tang X-H, Roelofs WL (1987) Sex pheromone blend discrimination by male             lobes. Proc R Soc Lond B Biol Sci 213:249–277.
    moths from E and Z strains of European corn borer. J Chem Ecol 13:143–151.                30. Lei H, Riffell JA, Gage SL, Hildebrand JG (2009) Contrast enhancement of stimulus in-
 5. Linn CE, Young MS, Gendle M, Glover TJ, Roelofs WL (1997) Sex pheromone blend                 termittency in a primary olfactory network and its behavioral significance. J Biol 8:21.
    discrimination in two races and hybrids of the European corn borer moth, Ostrinia         31. Van Rullen R, Thorpe SJ (2001) Rate coding versus temporal order coding: What the
    nubilalis. Physiol Entomol 22:212–223.                                                        retinal ganglion cells tell the visual cortex. Neural Comput 13:1255–1283.
 6. Willis MA, Baker TC (1988) Effects of varying sex pheromone component ratios on the       32. Gawne TJ, Kjaer TW, Richmond BJ (1996) Latency: Another potential code for feature
    zigzagging flight movements of the oriental fruit moth, Grapholita molesta. J Insect           binding in striate cortex. J Neurophysiol 76:1356–1360.
    Behav 1:357–371.                                                                          33. Gollisch T, Meister M (2008) Rapid neural coding in the retina with relative spike
 7. Liu Y-B, Haynes KF (1992) Filamentous nature of pheromone plumes protects integrity           latencies. Science 319:1108–1111.
    of signal from background chemical noise in cabbage looper moth, Trichoplusia ni. J       34. Johansson RS, Birznieks I (2004) First spikes in ensembles of human tactile afferents
    Chem Ecol 18:299–307.                                                                         code complex spatial fingertip events. Nat Neurosci 7:170–177.
 8. Jones C (1983) On the structure of instantaneous plumes in the atmosphere. J Hazard       35. Hopfield JJ (1995) Pattern recognition computation using action potential timing for
    Mater 7:87–112.                                                                               stimulus representation. Nature 376:33–36.
 9. Riffell JA, Abrell L, Hildebrand JG (2008) Physical processes and real-time chemical      36. Margrie TW, Schaefer AT (2003) Theta oscillation coupled spike latencies yield com-
    measurement of the insect olfactory environment. J Chem Ecol 34:837–853.                      putational vigour in a mammalian sensory system. J Physiol 546:363–374.
10. Laurent G, Davidowitz H (1994) Encoding of olfactory information with oscillating         37. Junek S, Kludt E, Wolf F, Schild D (2010) Olfactory coding with patterns of response
    neural assemblies. Science 265:1872–1875.                                                     latencies. Neuron 67:872–884.
11. Laurent G, Wehr M, Davidowitz H (1996) Temporal representations of odors in an            38. Cleland TA, Johnson BA, Leon M, Linster C (2007) Relational representation in the
    olfactory network. J Neurosci 16:3837–3847.                                                   olfactory system. Proc Natl Acad Sci USA 104:1953–1958.
12. Laurent G (1999) A systems perspective on early olfactory coding. Science 286:723–728.    39. Zavada A, Buckley CL, Martinez D, Rospars J-P, Nowotny T (2011) Competition-based

13. Vickers NJ, Christensen TA, Baker TC, Hildebrand JG (2001) Odour-plume dynamics               model of pheromone component ratio detection in the moth. PLoS One 6:e16308.
    influence the brain’s olfactory code. Nature 410:466–470.                                  40. Linster C, Kerszberg M, Masson C (1994) How neurons may compute: The case of
14. Baker TC, Hansson BS, Löfstedt C, Löfqvist J (1988) Adaptation of antennal neurons in         insect sexual pheromone discrimination. J Comput Neurosci 1:231–238.
    moths is associated with cessation of pheromone-mediated upwind flight. Proc Natl          41. Kwok YC (2007) Encoding of odour blends in the moth antennal lobe. PhD thesis
    Acad Sci USA 85:9826–9830.                                                                    (University of Leicester, Leicester, UK).
15. Bhandawat V, Maimon G, Dickinson MH, Wilson RI (2010) Olfactory modulation of             42. Chen HT, Ng KT, Bermak A, Law MK, Martinez D (2011) Spike latency coding in a bi-
    flight in Drosophila is sensitive, selective and rapid. J Exp Biol 213:3625–3635.              ologically inspired microelectronic nose. IEEE Trans Biomed Circuits Syst 5:160–168.
16. Tumlinson JH, et al. (1989) Identification of a pheromone blend attractive to Man-         43. Christensen TA, Lei H, Hildebrand JG (2003) Coordination of central odor repre-
    duca sexta (L.) males in a wind tunnel. Arch Insect Biochem Physiol 10:255–271.               sentations through transient, non-oscillatory synchronization of glomerular output
17. Kaissling KE, Hildebrand JG, Tumlinson JH (1989) Pheromone receptor cells in the              neurons. Proc Natl Acad Sci USA 100:11076–11081.
    male moth Manduca sexta. Arch Insect Biochem Physiol 10:273–279.                          44. Ito I, Bazhenov M, Ong RC, Raman B, Stopfer M (2009) Frequency transitions in odor-
18. Christensen TA, Hildebrand JG (1987) Male-specific, sex pheromone-selective pro-               evoked neural oscillations. Neuron 64:692–706.
    jection neurons in the antennal lobes of the moth Manduca sexta. J Comp Physiol A         45. Börgers C, Kopell N (2003) Synchronization in networks of excitatory and inhibitory
    Neuroethol Sens Neural Behav Physiol 160:553–569.                                             neurons with sparse, random connectivity. Neural Comput 15:509–538.
19. Kanzaki R, Soo K, Seki Y, Wada S (2003) Projections to higher olfactory centers from      46. Kopell N, Ermentrout B (2004) Chemical and electrical synapses perform comple-
    subdivisions of the antennal lobe macroglomerular complex of the male silkmoth.               mentary roles in the synchronization of interneuronal networks. Proc Natl Acad Sci
    Chem Senses 28:113–130.                                                                       USA 101:15482–15487.
20. Haupt SS, Sakurai T, Namiki S, Kazawa T, Kanzaki R (2010) Olfactory information           47. Martinez D (2005) Oscillatory synchronization requires precise and balanced feedback
    processing in moths. The Neurobiology of Olfaction, ed Menini A (CRC, Boca Raton,             inhibition in a model of the insect antennal lobe. Neural Comput 17:2548–2570.
    FL), pp 71–112.                                                                           48. Martinez D, Montejo N (2008) A model of stimulus-specific neural assemblies in the
21. Lei H, Christensen TA, Hildebrand JG (2002) Local inhibition modulates odor-evoked            insect antennal lobe. PLoS Comput Biol 4:e1000139.
    synchronization of glomerulus-specific output neurons. Nat Neurosci 5:557–565.             49. Krofczik S, Menzel R, Nawrot MP (2009) Rapid odor processing in the honeybee an-
22. Waldrop B, Christensen TA, Hildebrand JG (1987) GABA-mediated synaptic inhibition             tennal lobe network. Front Comput Neurosci 2:9.
    of projection neurons in the antennal lobes of the sphinx moth, Manduca sexta. J          50. Linster C, Masson C, Kerszberg M, Personnaz L, Dreyfus G (1993) Computational diversity
    Comp Physiol A Neuroethol Sens Neural Behav Physiol 161:23–32.                                in a formal model of the insect olfactory macroglomerulus. Neural Comput 5:228–241.
23. Zhang HG, Lee HJ, Rocheleau T, ffrench-Constant RH, Jackson MB (1995) Subunit             51. Sivan E, Kopell N (2006) Oscillations and slow patterning in the antennal lobe. J
    composition determines picrotoxin and bicuculline sensitivity of Drosophila gamma-            Comput Neurosci 20:85–96.
    aminobutyric acid receptors. Mol Pharmacol 48:835–840.                                    52. Peron S, Gabbiani F (2009) Spike frequency adaptation mediates looming stimulus
24. Christensen TA, Waldrop BR, Hildebrand JG (1998) Multitasking in the olfactory sys-           selectivity in a collision-detecting neuron. Nat Neurosci 12:318–326.
    tem: Context-dependent responses to odors reveal dual GABA-regulated coding               53. Homberg U, Christensen TA, Hildebrand JG (1989) Structure and function of the
    mechanisms in single olfactory projection neurons. J Neurosci 18:5999–6008.                   deutocerebrum in insects. Annu Rev Entomol 34:477–501.
25. Khawaled R, Bruening-Wright A, Adelman JP, Maylie J (1999) Bicuculline block of           54. Skinner FK, Chung JYJ, Ncube I, Murray PA, Campbell SA (2005) Using heterogeneity
    small-conductance calcium-activated potassium channels. Pflugers Arch 438:314–321.             to predict inhibitory network model characteristics. J Neurophysiol 93:1898–1907.
26. Debarbieux F, Brunton J, Charpak S (1998) Effect of bicuculline on thalamic activity: A   55. Burges CJ (1998) A tutorial on support vector machines for pattern recognition.
    direct blockade of IAHP in reticularis neurons. J Neurophysiol 79:2911–2918.                  Data Mining and Knowledge Discovery 2:121–167.

Belmabrouk et al.                                                                                                    PNAS | December 6, 2011 | vol. 108 | no. 49 | 19795
Supporting Information
Belmabrouk et al. 10.1073/pnas.1112367108
SI Text                                                                                                        gγ (ms/cm2)                Vγ (mV)
Quantitative Coding. In the projection neuron (PN) model, the
stimulation current Istim reflects the excitatory drive from ol-       INa = gNAm h(V − VNA)
                                                                                                                 190                         50
                                                                      IK = gKn4(V − VK)                           60                       −100
factory receptor neuron (ORN) inputs (SI Text, Model). The PN
                                                                      ICa = gCas∞(V − VCa)                         5                        120
absolute latency computed as the time to first spike relative to
                                                                      IAHP = gAHPq(V − VAHP)                       4                       −140
stimulus onset was found to be inversely proportional to Istim        IL = gL(V − VL)                              0.1                      −67
(Fig. S4) (S1):

                                        1                               The parameters m, n, and h are, respectively, described by
                                 L∝         :                 [S1]    (Eqs. S7–S9)
                                                                                                  ¼ αm ð1 − mÞ − βm m;                      [S7]
From the simulations with varying concentration K and pro-                                     dt
portion P (Fig. 3 A and B), Istim ðBALÞ ¼ KPI0PN , Istim ðC15Þ ¼
Kð1 − PÞI0PN , and Istim ðgeneralistsÞ ¼ KI0PN , with I0PN ¼ 13:12                           dn
                                                                                                ¼ αn ð1 − nÞ − βn n ; and                   [S8]
μA=cm2 (SI Text, Model). Therefore, the absolute latency of the                              dt
three groups of PNs is well-described by the following functions
(Figs. S5 and S6) (S2–S4):                                                                     dn
                                                                                                  ¼ αh ð1 − hÞ − βh h:                      [S9]
                               LBAL ∝    ;                    [S2]    αγ and βγ (γ = m, n, and h) are listed below.

                         LC15 ∝            ; and              [S3]                               αγ                                  βγ
                                  Kð1 − PÞ
                                                                                          0:32ðV þ 54Þ                         0:28ðV þ 27Þ
                                             1                                                   V þ 54                           V þ 27
                             Lgeneralists   ∝ :               [S4]                     1 − exp −                             exp           −1
                                             K                                                      4
                                                                                         0:032ðV þ 52Þ                                 V þ 57
Expression S4 indicates that the latency of the generalist PNs is     n                                                    0:5exp −
                                                                                                 V þ 52                                  40
invariant to the proportion P while being sensitive to the con-                        1 − exp −
                                                                                                    5                                 4
centration K (behaving as 1/K).                                                                    V þ 50                               V þ 27
                                                                      h                0:128exp −                            1 þ exp −
                                                                                                      18                                  5
Qualitative Coding. By using expressions S2 and S3, the normal-
ized difference latency is expressed as (Eq. S5)

                               1     1
                                  −                                       The activation function for ICA is (Eq. S10)
              LC15 − LBAL Kð1 − PÞ KP
                          ∝             ¼ 2P − 1:             [S5]
              LC15 þ LBAL    1    1
                               þ                                                                              1
                            KP Kð1 − PÞ                                                      S∞ ¼                    :                   [S10]
                                                                                                                V þ 25
                                                                                                      1 þ exp −
Thus, the normalized difference latency is sensitive to the                                                       5
proportion P while being invariant to the concentration K, which
is in agreement with the simulations results (red curves in Fig. 3    Finally, the calcium-dependent activation function for IAHP is
A and B).                                                             described by (Eq. S11)

Nonoscillatory vs. Oscillatory Synchronization. Our macroglomerular                                          ½CaŠ
                                                                                                      q¼                                   [S11]
complex (MGC) model exhibited nonoscillatory synchronization,                                              30 þ ½CaŠ
which persisted when weak excitatory connections were added
                                                                      with (Eq. S12)
from PNs to local neurons (LNs) (Fig. S7, black and blue traces).
On the contrary, strong PN-LN excitatory connections elicited fast                &
                                                                          d½CaŠ      − 0:002ICA − 0:0000125½CaŠ for LN-lla
∼40-Hz oscillations (Fig. S7, red trace).                                       ¼                                                          [S12]
                                                                            dt     − 0:002ICA − 0:000125½CaŠ for LN-llb:
Model. LN model. The LN model is a conductance-based model with
intrinsic currents modified from ref. 1. The voltage V is described
by (Eq. S6)                                                           PN model. The projection neuron model is inspired by two models
                                                                      (2, 3), which we combined into one consistent model. The
                  dV                                                  membrane potential V is governed by (Eq. S13)
              C      ¼ Istim − INa − IK − ICa − IAHP − IL ;   [S6]
                  dt                                                           dV
                                                                           C      ¼ IStim − INa − IK − ICa − Isk − IA − IL − IGABA :       [S13]
where C = 1 μF/cm2 and the parameters are defined as below.                     dt

Belmabrouk et al.                                                                                 1 of 6
The ionic currents are described as follows (σ = 1).                                                                                      À                 Á
                                                                                                                  Istim2 ¼ Istim1 ðt0 Þexp − ðt − t0 Þ=τfall ;                 [S19]

                                                     gγ (mS/cm2)              Vγ (mV)
                                                                                            where t0 = 50 ms, τrise = 5 ms, and τfall = 4 s for all simulations
                                                                                            except for the repetition task (Fig. 2D), where τfall = 400 ms to
INa ¼ gNa m3 ð1 − WÞðV − VNa Þ
             ∞                                         120                       55         avoid overly long stimulations. The Imax value depends on the
IK = gK (W/σ)4 (V − VK)                                  5                      −70         neuron type and simulation task.
ICa = gCaX (V − VCa)                                     0.1                    124           For all simulations, except for the test of varying concentration
Isk ¼ gsk q2 ðV − VK Þ
           ∞                                             2                      −70         and component proportions (Figs. 3 and S3), we used the Imax
IA = gAA∞B (V − VK)                                     22                      −70         values detailed in the table below.
IL = gL(V − VL)                                          0.3                    −50
                                                                                                                                        Imax (for                  Imax (for
  The functions W, B, and X are defined by first-order differ-                                Type of                Imax (for         uniglomerular              multiglomerular
ential equations (here, γ = W, B, and X) (Eq. S14):                                         input                 LN; μA/cm2)         PN; μA/cm2)                PN; μA/cm2)

                                                                                            Single                     15                     35                      35
                                     dγ γ ∞ − γ
                                        ¼       ;                                [S14]        compound
                                     dt    τγ                                               Blend                      30                     35                      70

and for W, X, B, m, and A (Eq. S15),
                                                                                               The blend is expressed as Blend = K[P·bombykal (BAL) +
                                           1                                                (1 − P)·E,Z-11,13-pentadecadienal (C15)]. In case of variations
                         γ∞   ¼                       :                       [S15]
                                1 þ exp − 2aγ V − Vγ1=2                                     in concentration (Fig. 3A), we fix P to 0.66 and vary K from 1 to
                                                                                            4. For varying proportions (Fig. 3B), we fix K to 2 and vary P
The parameter for these functions is λ = 0.05. The calcium-                                 from 0.3 to 0.7 (I0PN ¼ 13:12 μA=cm2 and I0LN ¼ 5:625 μA=cm2 ).
dependent activation function for Isk is given by (Eq. S16)                                 In practice, Imax was changed according to the table below.

     q∞ ðCsk Þ ¼                                                    !;         [S16]
                                                            Csk − Cr                        Type of       Imax (for         Imax (for        Imax (for                Imax (for
                   1 þ exp − 1:120 − 2:508log                                               input           LN)             BAL-PN)          C15-PN)            multiglomerular PN)
                                                                                            Blend         KI0LN             KPI0PN          Kð1 − PÞI0PN               KI0PN
where Csk denotes the concentration in the sk calcium domain
and Cr (113 nM) is the resting calcium concentration.
                                                                                            GABAergic current IGABA. The inhibitory synaptic current, IGABA,
                Vγ 1=2                                                                      which is added to the current balance equation of the post-
       aγ       (mV)                                    τγ (ms)                             synaptic PN (4), is given by (Eqs. S20 and S21)

                                                                1                                                  IGABA ¼ gGABA sðV − VGABA Þ and                             [S20]
W    0.055       −46          τw ðVÞ ¼
                                         λexpðaw ðV − VW1=2 ÞÞ þ λexpð − aw ðV − Vw1=2 ÞÞ

B    −0.1       −70                                       20                                                        ds     À    Á
                                                                                                                       ¼ αT Vpre ð1 − sÞ − s=τGABA ;                           [S21]
X    0.071      −10                                       50                                                        dt
A    0.02       −20                                       —
m    0.065      −31                                       —
                                                                                            where s is the fraction of open channels. The release of neuro-
                                                                                            transmitter T was modeled by (Eq. S22)
                                                                                                                      À    Á                 1
    The calcium pools evolve according to (Eq. S17)                                                                  T Vpre ¼               À         Á;                       [S22]
                                                                                                                                     1 þ exp − Vpre =2
                         dCsk               1
                              ¼ − αsk ICa − ðCsk − Cr Þ;                         [S17]      with synaptic parameters VGABA = −75 mV, α = 6.25 ms−1,
                          dt               τsk
                                                                                            τGABA = 10 ms, gGABA = 0.08 mS/cm2 for LN-IIa, and gGABA =
with time constant τsk = 250 ms and Faraday scaling factor                                  0.04 mS/cm2 for LN-IIb.
αsk = 0.9.                                                                                  Kullback–Leibler distance. The Kullback–Leibler distance (KLD;
Stimulation input current IStim. In our model, the ORN output is                            also called the relative entropy) between two distributions p and
replaced by a current injection with the same experimentally                                q of a discrete random variable is expressed as (Eq. S23)
observed time shape. Each neuron in our MGC model will receive
an initial input when an odor is detected according to (Eq. S18)                                                                     X                 pðiÞ
                                                                                                                         KLD ¼              pðiÞ log        :                  [S23]
                         Istim1 ¼ Imax ð1 − expð − t=τrice ÞÞ:                   [S18]
                                                                                            In Fig. 1C, the KLD is computed between the real proportions p
After the odor was removed at t = t0, the neuron’s input is
                                                                                            and the proportions q obtained in simulations.
governed by (Eq. S19)

1. Benda J, Herz AV (2003) A universal model for spike-frequency adaptation. Neural         3. Av-Ron E (1994) The role of a transient potassium current in a bursting neuron model. J
   Comput 15:2523e2564.                                                                        Math Biol 33:71e87.
2. Roper P, Callaway J, Shevchenko T, Teruyama R, Armstrong W (2003) AHP’s, HAP’s and       4. Skinner FK, Chung JY, Ncube I, Murray PA, Campbell SA (2004) Using heterogeneity
   DAP’s: How potassium currents regulate the excitability of rat supraoptic neurons. J        to predict inhibitory network model characteristics. J Neurophysiol 93:1898e
   Comput Neurosci 15:367e389.                                                                 1907.

Belmabrouk et al.                                                                                                                    2 of 6
Fig. S1. Schematic of the MGC network model having 86 LNs-IIa, 68 LNs-IIb, and 41 PNs. Both LN types are multiglomerular, and 30% of the PNs are mul-
tiglobular. The remaining 70% are uniglomerular, connected to either BAL or C15 ORNs (cumulus or toroid glomerulus) with equal probability. Inhibitory LNs
are randomly connected to the PNs with probabilities p(LN-IIa → PN) and p(LN-IIb → PN). Unless stated otherwise, p(LN-IIa → PN) = 0.5 and p(LN-IIb → PN) = 0.2.
Note that the model is essentially feed-forward, lacking excitatory connections from PNs to LNs. Adding excitatory connections from PNs to LNs, however, did
not dramatically change the results (Discussion).

Belmabrouk et al.                                                                                              3 of 6
                                  A                                               B

                                  D                                               C

Fig. S2. Effects of removing LN-IIa or LN-IIb. (A) Without LN-IIa, the proportions of response patterns were obtained as in Fig. 1 (30 network simulations of 1 s
biological time). The biphasic patterns represented 68% of the responses, which corresponds to the proportions of biphasic plus triphasic patterns in the intact
network (Fig. 1). LNs-IIa are, therefore, responsible for the first inhibitory phase I1 in the triphasic response. (B) Without LN-IIb, the duration of the second
inhibitory phase I2 reduced significantly: 106 ± 2.68 ms without LN-IIb vs. 112 ± 5 ms for the intact network (P < 0.05, Kruskal–Wallis test). The LNs-IIb are,
therefore, involved but not alone responsible for I2. (C) Blend-enhanced effect on synchrony in the intact MGC network without LN-IIa or LN-IIb. Synchrony was
expressed as the percentage of synchronous spikes found in all pairs of specific PNs (Methods). In all cases (intact MGC or without LN-IIa or LN-IIb), synchrony
was higher with the blend than with single compounds (P < 0.05, Kruskal–Wallis test). Nevertheless, the blend-enhanced effect on synchrony was more
pronounced in the intact network (69 ± 27% with the blend vs. 59 ± 26% with single compounds) than without LN-IIa or LN-IIb. Although both LN types are
involved in PN synchrony, removing the LNs-IIa had significantly more impact (69 ± 27% for the intact MGC vs. 51 ± 23% without LN-IIa and 56 ± 28% without
LN-IIb, P < 0.05, Kruskal–Wallis test followed by Mann–Whitney test). (D) Representative voltage traces of two PNs (red and black traces) in the intact network
and without LN-IIa. In the intact network (Upper), inhibition I1 before E plays the role of a reset by relaxing the PN activities to a hyperpolarized state. As
a consequence, the PNs fired synchronously when inhibition vanished. In the network without LN-IIa (Lower), I1 was not strong enough to play the role of
a reset and eliminate the influence of initial conditions.

Belmabrouk et al.                                                                                                4 of 6

                    B                                         C                                     D

Fig. S3. Responses of the MGC model to blend-like stimuli in different proportions P and concentrations K. The binary mixture model is K(P × BAL + (1 − P) ×
C15). (A) With P = 2/3 (proportion of BAL two times as big as C15), the spike raster activities revealed three distinct populations of synchronized neurons. The
generalist PNs (black) fired first as they received inputs from both glomeruli and were followed by the BAL-PNs (blue) and C15-PNs (red). (B) Mean absolute
latencies relative to the stimulus onset were computed for the three populations of neurons in each trial and plotted in a 3D space. The points in red cor-
respond to the correct proportion (80 samples for proportion P = 2/3 and concentration K = 1, 2, 3, 4). The points in blue correspond to incorrect mixture ratios
(100 samples for P = 0.3, 0.4, 0.5, 0.6, 0.7 and K = 1). We trained a linear support vector machine (SVM) to discriminate the correct proportion P = 2/3 against
incorrect mixture ratios (red vs. blue points). The mean discriminant plane depicted in the figure is obtained by cross-validation (resampling 10 times the
training and test sets and averaging the results; discussed in Methods). (C) Same procedure as in B but for the mean firing rates estimated in a window of 200
ms after stimulus onset. (D) Classification performances were higher for latency than for rate (99.7 ± 0.5% of success with latency vs. 96.4 ± 1.3% with rate,
P < 0.05, Kruskal–Wallis test).

Fig. S4. Absolute latency L (ms) of the three groups of PNs vs. stimulus intensity Istim (μA/cm2). Values for Istim were computed from the simulations with
varying concentrations K and proportions P (Fig. 3 A and B) as Istim ðBALÞ ¼ KPI0PN , Istim ðC15Þ ¼ Kð1 − PÞI0PN , and Istim ðgeneralistsÞ ¼ KI0PN , with
I0PN ¼ 13:12 μA=cm2 (SI Text). Mean and SDs for latencies were estimated over 20 runs. The black, blue, and red points are for the generalist PNs, C15-PNs, and
BAL-PNs, respectively. The dashed curves represent the fits ∝1=Istim : LBAL;C15 ¼ 706=Istim þ 38ms and Lgeneralists ¼ 152=Istim þ 39ms.

Belmabrouk et al.                                                                                                5 of 6
Fig. S5. Absolute latency L (ms) of the three groups of PNs vs. stimulus concentration K (proportion P = 2/3). Mean and SDs were estimated over 20 runs. The
dashed curves represent the equations LC15 ¼ 706=ðKð1=3ÞI0PN Þ þ 38ms, LBAL ¼ 706=ðKð2=3ÞI0PN Þ þ 38ms, and Lgeneralists ¼ 152=ðKI0PN Þ þ 39ms, with
I0PN ¼ 13:12 μA=cm2 (Fig. S4).

Fig. S6. Absolute latency L (ms) of the three groups of PNs vs. stimulus proportion P (concentration K = 2). Mean and SDs were estimated over 20 runs. The
dashed curves represent the equations LBAL;C15 ¼ 706=ð2PI0PN Þ þ 38ms, LC15 ¼ 706=ð2ð1 − PÞI0PN Þ þ 38ms, and Lgeneralists ¼ 152=ð2I0PN Þ þ 39ms, with
I0PN ¼ 13:12 μA=cm2 (Fig. S4).

Fig. S7. Power spectra of averaged PN activity for different excitatory synaptic conductance gPN→LN from PNs to LNs (probability of connection = 0.8). The
excitatory synaptic current was modeled as a sum of exponential decays gPN→LN i;f e − ðt − ti Þ=τ , where tif is firing time f of the ith presynaptic LN and τ ¼ 10ms.

The power spectrum had a sharp peak near 40 Hz for strong excitatory connections (gPN→LN ¼ 15, red trace).

Belmabrouk et al.                                                                                                    6 of 6

Shared By: