## Abstract

The spatial and temporal characteristics of the visual and acoustic sensory input are indispensable attributes for animals to perform scene analysis. In contrast, research in olfaction has focused almost exclusively on how the nervous system analyzes the quality and quantity of the sensory signal and largely ignored the spatiotemporal dimension especially in longer time scales. Yet, detailed analyses of the turbulent, intermittent structure of water- and air-borne odor plumes strongly suggest that spatio-temporal information in longer time scales can provide major cues for olfactory scene analysis for animals. We show that a bursting subset of primary olfactory receptor neurons (bORNs) in lobster has the unexpected capacity to encode the temporal properties of intermittent odor signals. Each bORN is tuned to a specific range of stimulus intervals, and collectively bORNs can instantaneously encode a wide spectrum of intermittencies. Our theory argues for the existence of a novel peripheral mechanism for encoding the temporal pattern of odor that potentially serves as a neural substrate for olfactory scene analysis.

- bursting olfactory receptor neuron
- interval estimation
- neural coding
- olfactory scene analysis
- point process
- uncoupled oscillators

## Introduction

Vision and audition provide dynamic perception of the world, commonly referred to as visual and auditory scene analysis, respectively, with enormous advantage for survival in higher vertebrates (Bregman, 1994). Given that olfaction is phylogenetically the oldest sensory system (Ache, 1991; Hildebrand, 1995), it would be surprising if animals that rely primarily on chemosensory information to quantify their sensory worlds did not evolve some version of olfactory scene analysis as an edge for survival (Mountain and Hubbard, 2001). Just as binocular disparity and interaural time differences provide critical information about the spatial perception in vision and audition, respectively, the irregular interval of odor encounters inherent in the odor plume structure, the intermittency, is a strong candidate for that purpose in olfaction (Moore and Atema, 1988; Murlis et al., 1992; Gomez et al., 1999).

Most research on olfactory coding has been devoted to the detection and discrimination of odor quality (Hopfield, 1999; Bazhenov et al., 2001; Bathellier et al., 2008; Raman et al., 2010). The study of intermittency in olfaction primarily has focused on the rapid discontinuous stimulus acquisition process [e.g., sniffing in mammals (Wachowiak, 2011), flicking in crustaceans (Schmitt and Ache, 1979)], or rapid stimulus flux (Brown et al., 2005; Geffen et al., 2009) that occurs over fractions of a second, in contrast to the second to many tens of seconds intermittency that characterizes the natural odor signal profile (Vickers et al., 2001; Webster and Weissburg, 2001; Gardiner and Atema, 2010). Because intermittency is a major cue in both passive and active sensing of the olfactory environment, one can hypothesize that there is a specialized subsystem that represents some form of the timing history of odor encounters.

Internal representation of interval timing is usually associated with higher-order brain function. Various central neural mechanisms, including pacemaker-accumulators, neural oscillators, and network dynamics have been proposed (Miall, 1989; Buhusi and Meck, 2005; Laje and Buonomano, 2013). Peripheral strategies, such as a system of uncoupled oscillating detectors, generally have not been considered. However, the potential for such a peripheral sensory system capable of dynamic information extraction exists in the olfactory system. Although the majority of olfactory receptor neurons (ORNs) in animals are canonical tonically active cells that respond to odor with changes in discharge rate, some are intrinsically or conditionally rhythmically active ORNs (Sicard, 1986; Frings and Lindemann, 1988; Reisert and Matthews, 2001a,b; Arnson and Holy, 2011) referred to here as “bursting” ORNs (bORNs). bORNs have been well characterized in the lobster olfactory organ, where they generate intrinsic bursts in response to brief odor stimuli in a phase-dependent manner (Bobkov and Ache, 2007). Although bORNs may potentially provide a basis for detecting odor intermittency, there are several fundamental questions that need to be addressed to ultimately understand whether bORNs actually have the potential to serve this purpose.

In the present study, we combine electrophysiological and theoretical modeling approaches to explore the tuning properties of individual lobster bORNs, and show that a heterogeneous population of such bORNs can encode a wide spectrum of stimulus intermittencies. Our findings argue for a heretofore unsuspected role for ORNs; encoding of the temporal structure of odor signals in a manner that provides a potential neural substrate for olfactory scene analysis.

## Materials and Methods

##### Preparation, electrophysiology, calcium imaging, experimental data analysis.

Experiments were conducted using a lobster ORN *in situ* preparation developed earlier (Bobkov and Ache, 2005; Ukhanov et al., 2011). A single annulus was excised from the olfactory organ (the lateral antennular filament) of lobsters of either sex, and the cuticle on the side opposite the olfactory (aesthetasc) sensilla was removed to provide better access to the cell bodies of the ORNs. After enzymatic (∼1 mg/ml trypsin, papain, collagenase for 10 min) treatment and mechanical clearing individual ORN somata were available for electrophysiological recording. The specimens were mounted on a plastic 35 mm Petri dish and placed on the stage of inverted microscopes (Axiovert 100, Zeiss). The cell bodies of the ORNs were continuously superfused with *Panulirus* saline (PS) containing the following (in mm): 486 NaCl, 5–14 KCl, 13.6 CaCl_{2}, 9.8 MgCl_{2} and 10 HEPES, pH 7.9. A second superfusion flow of PS allowed odorants to be delivered exclusively to the olfactory sensilla containing the outer dendrites of the ORNs. Both superfusion contours were gravity fed at constant rates of flow. The odorant was an aqueous extract of Tetra-Marine (TET), a commercially available marine fish food. Flakes of TET were powdered, dissolved in water (0.1 g/ml), filtered through a 0.2 μm syringe filter, and diluted 1:200 in PS for experiments. The total maximum concentration applied was estimated to be ∼60–100 μm. The concentration of any individual component might be expected to be 1 nm-100 nm range, depending on total number of components, consistent with concentrations able to elicit chemotaxis in spiny lobsters (Ache et al., 1978). The odorant stream was switched with the flow of PS that otherwise continuously superfused the sensilla (both ∼250 μl/min) using a multichannel rapid solution changer (RSC-160, Bio-Logic).

Electrophysiological evaluation of ORN activity was performed using both extracellular loose-patch and whole-cell recording. Patch electrodes were pulled from borosilicate capillary glass (Sutter Instrument BF150-86-10) using a Flaming-Brown micropipette puller (model P87). Resistance of the electrodes was 1–5 mΩ as measured in PS. For whole-cell recordings, the electrode contained intracellular solution: KCl 180 or 210 mm, NaCl 30 or 0 mm, GTP 0.5 mm, ATP 0.5 mm, MgCl_{2} 1 mm, Glucose 696 mm, HEPES 10 mm, and Tris-base to adjust pH to 7.8.

For calcium signal-based evaluation of ORN activity, the antennular segments were placed in an Eppendorf tube in PS containing the fluorescent calcium indicator (Fluo-4AM) 10 μm prepared with 0.2–0.06% Pluronic F-127 (Invitrogen). The tube was shaken for 30–60 min on an orbital shaker (∼70 rpm). The tissue was transferred into fresh PS and mounted for imaging. Fluorescence imaging was performed on an inverted microscope (Olympus IX-71) equipped with a cooled CCD camera (ORCA R2, Hamamatsu) under the control of Imaging Workbench 6 software (INDEC Systems). The stimulation paradigm was similar to that used in electrophysiological experiments. A standard FITC filter set (excitation at 510 nm, emission at 530 nm) was used. Images were collected at the rate of ∼4.23 Hz. Recorded data were stored as image stacks, and analyzed off-line using Imaging Workbench 6 or ImageJ 1.42 and pCLAMP.

Voltages/currents were measured with an Axopatch 200B patch-clamp amplifier (Molecular Devices) using an AD–DA converter (Digidata 1320A, Molecular Devices), low-pass filtered at 5 kHz, sampled at 5–20 kHz. Data were collected and analyzed with pCLAMP 9.2 software (Molecular Devices) in combination with Sigma Plot 10.0 (SPSS). The time of occurrence of the spike was taken as the time of peak current deflection, i.e., the peak of the spike. Burst detection and analysis was performed using the burst analysis protocol provided in pCLAMP 9.0. Burst delimiting interval and a minimum number of spikes in a burst were individually specified for each ORN. Interburst intervals (IBIs) were usually taken as the time between the first spikes of two subsequent bursts. The bin width is varied in the range 0.25–2 s depending on individual bORN bursting frequencies. The data are presented as the mean ± SEM of *n* observations unless otherwise noted. All recordings were performed at room temperature.

##### Probabilistic modeling.

We model the bORN bursts as a renewal process with external resets (see Fig. 3). Let *F* be the cumulative distribution of the spontaneous interburst intervals. Under the renewal assumption, the dependence of the past state quickly is washed out, enabling simple approximations. In other words, when the oscillator is noisy, no matter if and when it bursted recently, the probability of bursting at a distant future *t* approaches a constant. In fact, one can show that the asymptotic distribution of the phase ϕ, defined as the time since last burst at time *t*, does not depend on *t* and is given by:
where μ is the mean of *F* (Karlin and Taylor, 1975). The phase distribution *F _{t}* (ϕ) converges to

*F*

_{∞}(ϕ) geometrically fast. We parameterize F(ϕ) as a normal distribution with mean μ and variance σ

^{2}, truncated to the positive reals.

Odor stimulus evokes a burst with probability *P _{e}*(ϕ), where ϕ is the time interval from last burst. We parameterize the evoked probability as a sigmoid

*P*(ϕ) =

_{e}**Y**

_{1}be the Bernoulli random variable for the response to the first stimulation of the bORN in an asymptotic state (taking a value of 1 for response, 0 for no response). The baseline response of a bORN is the probability of response,

*p*= Pr[

**Y**

_{1}= 1]. By taking the expectation on the (asymptotic) time since last burst, we obtain: Note that the empirically estimated probability of the baseline has a positive bias of μΔ

*t*, if the actual window for burst detection is Δ

*t*, because the stimulation-induced response and spontaneous bursting are virtually identical (see Fig. 2

*b*).

##### Most recent interval tuning curve.

We now derive the response probability of a bORN in the asymptotic state exposed to a pair of stimulations. Let **Y**_{2}|Δ be the Bernoulli random variable for the response to the second stimulation Δ seconds after the first stimulation to the bORN. We assume the bORN was in the asymptotic state at the time of the first stimulus. If *q*(Δ) = Pr [**Y**_{2} = 1 | Δ] is different from *p*, it indicates that the system response is tuned to the stimulation interval Δ. Using the law of total probability conditioned on the first response, we have:
where Pr[**Y**_{1} = 0] is simply equal to 1 − Pr[**Y**_{1} = 1] = 1 − *p*. Note that Pr[**Y**_{2} = 1|Δ, **Y**_{1} = 0] ≠ Pr[**Y**_{1} = 1], because if the bORN did not respond to the first stimulus, it is not in the asymptotic state anymore; it selectively lost some portion of the phase distribution due to the first stimulus. Taking the expectation over the time since the most recent burst ϕ,
The distribution of the time since last burst ϕ depends on how many spontaneous bursts were there in the interval Δ. The distribution of *k*′th spontaneous burst given **Y**_{1} = 1 corresponds to the *k*-fold convolution of *F*, denoted as *F _{k}*. The phase distribution at Δ given

*k*spontaneous bursts, but not

*k*+ 1 is given by the product of

*k*-fold convolution and the probability of having no burst in Δ − ψ: If

**Y**

_{1}= 0, right after the first stimulation, the phase distribution is proportional to (1 −

*P*(Δ))

_{e}*F*

_{∞}(Δ). This quickly converges to the asymptotic distribution

*F*

_{∞}in the order of (μ −

*x*

_{0}) seconds. Therefore, we approximate Pr[

**Y**

_{2}= 1|δ,

**Y**

_{1}= 0] ≈

*p*. Finally, we obtain

*q*(Δ), the (approximate) most recent interval tuning curve; the conditional probability of firing given that the last odor stimulation was Δ seconds ago, as in Eq. 6.

##### Experimentally measured parameters and population extrapolation.

Each bORN is modeled with four parameters: (μ, σ) parameter of ϕ which corresponds to the mean and SD of the normal distribution before truncation, and (*x*_{0}, *b*) for shift and slope of the sigmoid *P _{e}*(ϕ).

We extrapolate a collection of sample parameter pairs [187 bORNs for μ, 28 for (μ, σ), 15 for (μ, σ, *x*_{0}, *b*)] obtained from electrophysiology to a large realistic population of bORNs each modeled with four parameters. Figure 5 shows a particular instantiation of the simulated parameters extrapolated from data as described below.

The mean parameter for the interburst interval (IBI) μ is strictly positive, and cannot be shorter than either the minimum burst duration nor the absolute refractory period of the bORN. Based on the observation from several bORNs, we consider this minimum to be 50 ms. The Weibull distribution (with 50 ms shift) fit the data best (Fig. 1*c*).

Experimentally, it is easier to detect bORNs with a shorter IBI. The procedure to find these bORNs was to patch random cells and wait to see if it bursts spontaneously or in an odor-induced manner. Hence, this distribution is likely biased toward shorter IBIs.

Once the mean IBI was determined, the SD of the IBI distribution was obtained by a linear fit. The variance was not constant, hence the fit was crosschecked with a generalized least-squares method, and the difference was negligible. The noise was assumed to be Gaussian with zero mean and SD estimated from the residuals. A sigmoid function with slope parameter 1/10 was used to model the change of variance (see Fig. 5*a*). The coefficient of variation, defined as the SD over the mean of the interval distribution, is ∼0.3. Hence, the point process is sub-Poisson as expected, and fairly regular.

The shift *x*_{0} in *P _{e}*(ϕ) was also linearly fitted (see Fig. 5

*b*). In this case, the bias was fixed to zero; a free bias results in negative

*x*

_{0}, which is unacceptable because it would imply that there is no refractory period.

The sigmoid slope parameter *b* in *P _{e}*(ϕ), is more tightly connected to

*x*

_{0}than μ (see Fig. 5

*c*). Because,

*x*

_{0}is generated from μ with noise, we assumed conditional independence of

*b*and μ given

*x*

_{0}. For small values of

*x*

_{0}, the slope of the regression line is flatter. In addition, the value is strictly positive; hence, we fit an exponential function. To generate strictly positive noise, zero mean Gaussian noise was added to the argument of the exponential. This induces slight positive bias to the generated points.

##### Calcium image data analysis.

Spontaneously oscillating cells were detected by inspecting average activity in the region of interests. Previous studies show that these oscillations correspond to bORN activity (Bobkov et al., 2012). The time series was denoised by removing the population average from each session (simultaneous recording of several cells). Then, low-frequency fluctuations were removed by low-pass filtering (cutoff 0.012 Hz) and subtracting. Three neurons were discarded from the analysis because calcium signal oscillations corresponding to bursting events weakened in amplitude over time making impossible to reliably detect bursts toward the end of the recording sessions. Burst detection was done by selecting a threshold crossing in the numerical-derivative. It was semiautomated: a Gaussian mixture model was fit to the distribution of the derivative to find the boundary, and then the threshold was manually fine-tuned. Evoked bursts were detected if there was a burst in the window after the stimulus with 1 s offset and size of 1–5 s.

Maximum a posteriori estimate of the parameters for the sigmoid *P _{e}*(ϕ) were obtained by imposing a Gaussian penalty of −λ‖1/

*b*‖

^{2}to the log-likelihood where λ = 0.01. Due to the limited number of trials (median of 7), estimated

*q*(Δ)'s were sometimes highly variable. Because the log-likelihood Eq. 8 is most sensitive to cells with extreme conditional probability, small errors in estimating

*q*(Δ) induce high variability in the decoding. To avoid the extreme variability, six cells with

*q*(Δ)'s <0.1 or >0.9 were excluded from the analysis.

## Results

As discussed earlier, intermittency estimation is likely a key component in olfactory scene analysis. Our aim is to show that a population of bORNs encode arbitrary odor intermittency (defined as the most recent interval between odor encounters), and moreover show that its instantaneous population activity can be linearly decoded. First, we describe the properties of bORNs via electrophysiology to show that (1) they are noisy oscillators with stereotypical burst patterns, and (2) they respond to odor in a phase-dependent manner. Based on the neurophysiology, we build a mathematical model that captures those two key observations. From the model, we found that the bursting response is tuned to the interval between odor stimulations. We demonstrate that a heterogeneous population of bORNs can flexibly encode a wide range of intermittencies, and propose a simple linear decoding scheme. Finally, we test the computational modeling by decoding intermittency from two sets of experimental data.

### Physiological characterization of bORN activity

As mentioned above, this study is based on a recent finding that the lobster olfactory organ contains at least two different subpopulations of ORNs (Bobkov and Ache, 2007; Ukhanov et al., 2011; Bobkov et al., 2012). Each subpopulation shows a different pattern of spontaneous activity and stimulus responsiveness. One subpopulation, heretofore unrecognized, shows spontaneous quasi-rhythmic bursts of action potentials (Fig. 1*a*). Using extracellular spike recording, we first characterized the spontaneous activity of rhythmically active ORNs (bORN). As shown for five bORNs, each cell has its own spontaneous rhythm of activity (burst frequency; Fig. 1*a*,*c*) and can be characterized by a specific burst pattern (Fig. 1*b*). The bursting parameters were largely consistent for any given cell over the time course of an experiment (Fig. 1*b*,*c*). Overall spontaneous burst frequency distribution (Fig. 1*d*) estimated from 187 recorded bORNs shows a wide range, well fit by a Weibull distribution.

### Odor-evoked activity of bORNs

To study the odor-evoked activity, the cells were stimulated by brief stimulus pulses of the same intensity (Fig. 2*a*; experimental trials for one bORN). The odor pulses trigger evoked bursts (Fig. 2*c*, middle) that are structurally similar to spontaneous bursts (Fig. 2*c*, top). Under whole-cell, voltage-clamp conditions (Fig. 2*d*), some bORNs were also able to generate both spontaneous (Fig. 2*d*, top) and odor-activated (Fig. 2*d*, middle) excitatory inward currents. As with the activity patterns recorded extracellularly, the burst-induced current patterns were similar (Fig. 2*d*, bottom).

We estimated the phase dependence of the evoked burst probability (Fig. 2*b*, blue symbols and line) as the ratio between the number of bursts occurring after stimulus presentation and the total number of stimulus presentations (usually 40–250 presentations per each ORN recorded extracellularly; Fig. 2*b*). For example, as shown for one bORN (Fig. 2*b*), the probability that the cell would generate a burst of action potentials upon stimulation ∼3 s after the last burst is 50% and almost 100% at a 6 s interval of quiescence. The spontaneous burst probability distributions (interburst interval histograms; Fig. 1*c*) can be expressed as cumulative distributions for comparison (Fig. 2*b*, red solid line). The interaction between the two sigmoidal curves determines the response characteristics of a specific bORN where an intermediate range of odor intermittency gives rise to the most reliable burst. Intuitively, a burst response is likely when stimulated during a longer quiescence, but at the same time, spontaneous bursts prevent such long inter burst intervals. Thus, discharge of the bursting cells is linked to the phase of bursting cycle relative to when the odor arrives, conferring the cell population with the potential to specifically encode the temporal structure of the odor stimulus.

### Computational modeling of bORNs

Because the functional hypotheses we address in this study involve the ensemble behavior of a population of bORNs rather than a single bORN activity, we used computational modeling and analysis to characterize it. Our primary goal is to derive the population-encoding model for temporal pattern of odor encounter and show that it is appropriately tuned to encode odor intermittency. Here we derive the conditional probability of the population-firing pattern given the temporal structure of the preceding odor encounter.

There are two major approaches to modeling pacemaker-like neuronal behavior: dynamical neuron modeling and point process modeling. Because the bORNs behave as stochastic oscillators and we are interested in the coding characteristics of a virtually independent ensemble (Bobkov et al., 2012), but not in the detailed mechanisms of bursting *per se*, a stochastic point process modeling approach is pursued. Specifically, we use a renewal process model where each burst generated by the bORN is considered as a unitary event, and the interval between the burst events is modeled as an independent draw from a distribution. Therefore each burst in the modeling illustrations (Fig. 2) is depicted as a symbol.

Let *F* be the common (cumulative) distribution from which the spontaneous interburst intervals are drawn for a specific bORN (Figs. 2*b*, 3*b*, red solid lines). If *F* is tightly concentrated around its mean, then the bORN is more regular (i.e., quasiperiodic), whereas a dispersed *F* will result in less predictable spontaneous interburst interval pattern. We model *F* as a truncated Gaussian distribution fitted from the data (Figs. 1*c*, 3*b*), because the time intervals are strictly positive. The other important quantity in our model is *P _{e}*(ϕ), the probability of odor-evoked bursting given that the last burst occurred ϕ seconds before the odor presentation. Due to refractoriness,

*P*(ϕ) is generally a monotonically increasing function and we parameterize its shape as a sigmoid function. Our bORN model for each cell is completely described by

_{e}*F*(·) and

*P*(·) (Fig. 3

_{e}*b*).

To model the stimulus-evoked burst generation, we modified the standard renewal model such that when stimulus is presented, with probability *P _{e}*(ϕ) it would burst and reset the phase. That is to say, the time interval in which the stimulus occurred becomes equal to the phase ϕ with probability

*P*(ϕ) (Figs. 2

_{e}*b*, 3

*b*, blue lines). From the sigmoidal shape of the conditional response probability

*P*(ϕ), one might expect that the neuron is more responsive to longer stimulus intervals. However, if the stimulus interval is long, a spontaneous burst is likely to happen in between, resulting in a small ϕ, and in turn, less responsiveness. Indeed, the interaction between

_{e}*P*(ϕ) and

_{e}*F*(ϕ) determines the tuning property of the bORN model to stimulus intermittency (Fig. 3

*c*).

### Theoretical tuning properties of bORNs

To investigate the change in response to odor depending on the intermittency, we derive the probability of response given a stimulus interval. The probability of response to a stimulus shortly after observing a burst Δ seconds ago is approximately the probability of having no spontaneous activity for Δ seconds and evoking a burst when the phase is Δ, that is, *P _{e}*(Δ) (1 −

*F*(Δ)). This approximation is only valid for small intervals because it ignores the possibility of having spontaneous bursts in the interval, and still responding to the stimulus. A better approximation is given by: where

*F*is the distribution of the

_{k}*k*′th burst which is a

*k*-fold convolution of

*F*, and

*c*). Note that lim

_{Δ→∞}

*q*(Δ) =

*p*, because the distribution of ϕ approaches

*F*

_{∞}.

From this analysis we can see that the bORN has a specific tuning characteristic for stimulus interval Δ. When *q*(Δ) is larger than *p*, the bORN response is amplified, and when it is smaller, the bORN response is reduced (Fig. 3*c*). The peaks of *q*(Δ) above *p* are analogous to natural resonance periods.

### Differential response to periodic stimulation

In natural environment, odor plumes produce a highly irregular temporal pattern of stimuli. However, to intuitively understand how bORNs respond to temporal patterns of stimulation, we demonstrated the response to an artificial periodic pattern of stimuli. A pair of subpopulations with identical spontaneous parameters, but with different *P _{e}*'s are set up; the difference in resulting tuning characteristic

*q*(Δ) can induce frequency selective partial synchronization (Fig. 4). The two homogeneous subpopulations in the absence of odor stimulus produce random and uncorrelated population activity. When an odor stimulus is presented, the population synchronize differentially and desynchronize over time due to their variability across intervals. Figure 4 shows subpopulation 1 amplifying for more sparse odor intermittency, whereas subpopulation 2 amplifying for the fast stimuli, which is also observable from the firing rate evolution over time. Note that due to the probabilistic nature of the response, even if all bORNs had the same period, they would not be perfectly entrained to the periodic stimulation. This clearly shows that the response of bORNs contains information regarding the history of input but in a probabilistic manner. In the next section, we show that a plausible heterogeneous population of bORNs reliably encodes the irregular temporal information of the stimulus.

### Decoding the stimulus interval from a heterogeneous population of bORNs

Heterogeneous populations of oscillators have been proposed to form a time basis in the brain (Miall, 1989; Matell and Meck, 2004). Miall (1989) showed that after full population synchrony, a subpopulation synchrony occurs when the least common multiple of periods is reached, and hence proposed beat (interference between frequencies that modulates amplitude) as a mechanism to form a set of oscillators with desired period from a smaller set of pacemaker neurons. However, this idea is not directly applicable to the population of bORNs, because two key assumptions are violated: (1) neurons are not perfect oscillators, and (2) full synchrony is not induced when the stimulus is given. Therefore, we extend this framework by relaxing these assumptions and proposing a probabilistic model for population coding of the temporal structure of the sensory input. From the population model, we will show that the instantaneous population response encodes sufficient information to infer the temporal structure of olfactory stimulus. Moreover, we also show that a simple maximum likelihood approach can reliably decode the stimulus time interval.

Although bORNs can be sharply tuned to stimulus intervals, each bORN has limited information content due to the uncertainty in both spontaneous and evoked bursting. Hence, a population is required to reliably encode the full spectrum of stimulus intervals. Suppose we have *N* independent bORNs each characterized by a unique pair of *F* and *P _{e}*. In addition, as can be seen from Figure 3

*c*, there are multiple stimulus intervals that correspond to the same probability of bursting. Therefore, decoding will be ambiguous for a homogeneous population; for unambiguous representation and decoding, a heterogeneous population is necessary.

Given the conditional probability of bursting *q _{i}*(Δ) for each bORN (indexed by

*i*), we define the population likelihood spectrum for decoding: where

*A*(and

*A′*) is the set of indices for the bORNs that responded to the current stimulation (or not). Therefore, estimation of an unknown Δ given the population response can be achieved by maximum likelihood (ML) estimation, that is, by finding the maximizing Δ in the log-likelihood function given by the following: where

*s*= 1, if neuron

_{i}*i*bursted, or

*s*= 0 otherwise, and

_{i}*c*(Δ) =

*q*(Δ)) does not depend on the response. We can design a neurally plausible decoding population that performs ML decoding where each neuron has a preferred intermittency Δ

_{i}*. Each decoding neuron receives a linear projection from the bORN population, and is laterally inhibited to implement a winner-take-all strategy. Because ML decoding first linearly projects the population response vector*

_{j}**= [**

*s**s*] to a set of weights

_{i}**(Δ**

*w**) = [log (*

_{j}*q*(Δ) − log (1 −

_{i}*q*(Δ))], and then takes the maximum of the projections adjusted by a bias

_{i}*c*(Δ

*).*

_{j}We tested ML decoding on a simulated population for a random stimulation timing pattern (see Fig. 6). To model the bORN ensemble behavior, we constructed an artificial population of bORNs based on experimentally obtained parameters (Fig. 5). The functional relation between the model parameters were extrapolated and used to construct an artificial population of bursting ORN models (see Materials and Methods). The population size was chosen to be *N* = 2000 based on a coarse estimate (∼2%) of the bORN population in a single lobster antennule.

Population-wise partial synchrony (which is not maintained due to heterogeneity) occurs as a response to each odor stimulus (Fig. 6*a*,*b*), and it can in turn be used to form the log-likelihood spectrum for decoding (Fig. 6*c*). The decoding recovers the stimulation interval with small error. Further analysis shows that the error systematically increases as the stimulation interval gets longer due to the broader peak in the likelihood function (Figs. 6*c*, 7). The results of the analysis suggest a coding scheme based on the ensemble of rhythmically active neurons that potentially provides a simple neural mechanism for instantaneous neuronal memory about the time since the last sensory event encountered by the ensemble.

### Decoding intermittency *in vivo*

To experimentally confirm the population coding, we reconstructed two virtual populations of cells recorded via either electrophysiology or calcium imaging.

In the electrophysiologically obtained population, we have aggregated data from different single-cell recordings to mimic a simultaneously observed population. Specifically, we probed 11 bORNs multiple times with a fixed timing pattern of two odor stimuli. From the burst time series, a realistic population response was constructed by considering individual trials as a simultaneous observation from different neurons. We used 10–40 sweeps of 1 min trials, consisting of stimulation at 11.0 s and 31.7 s resulting in 210 trials (Fig. 8*a*). We stack those 210 trials to mimic a population response from 210 bORNs. From the observation of the second stimulation, a window of 0.5–2 s was selected, and all the bursts in the window were considered as the response to the stimulus. We compute the population log-likelihood spectrum, and the ML decoding result yields 23.2 s, compared with the true value of 20.7 s (Fig. 8*b*). The mismatch is within the range of expectation because decoding from a simulated bORNs with the same conditions gave similar results (19.1 ± 6.7), with the larger confidence interval presumably reflecting limited heterogeneity and small population size.

In the calcium imaging-based population, we have identified oscillating cells and simultaneously recorded their activity up to six neurons at a time. Using a similar procedure of stacking single trials and aggregating over recording sessions, we constructed a virtual population of 400 bORNs from 59 original neurons (Fig. 9*a*,*b*). Although we assumed independence when aggregating different trials and sessions, within each session the cells are not necessarily independent. The ML decoding resulted in 21.2 s, compared with the true interval of 20 s (Fig. 9*c*).

## Discussion

Using electrophysiological, calcium imaging, and computational modeling approaches, we identified a functional neural model with the potential to capture signal intermittency. The proposed neural code is based on a rhythmically active subset of ORNs referred to as bORNs that have the interesting capacity to encode the temporal properties of intermittent odor signals. These neurons appear to function independently (Bobkov et al., 2012), with correlated activity induced by common sensory input. Based on the differences in the bORNs' inherent rate of bursting discharge and the phase dependency of their response to odor stimulation, each neuron selectivity modulates responses to a relatively narrow range of stimulus intervals.

Along with this finding, a more general hypothesis emerged from the computational analysis of bORN ensemble activity, which is that the instantaneous state of the heterogeneous population can reliably encode the time since the last stimulus. Although neuronal mechanisms for encoding and processing time intervals are still largely unknown, particularly in the seconds-to-minutes range as in the present study (Miall, 1989; Sumbre et al., 2008; Simen et al., 2011; Allman and Meck, 2012), the idea of using sensory modality-specific timing mechanisms instead of a centralized clock has been proposed for other systems (Miall, 1989; Buhusi and Meck, 2005; for review, see Bueti, 2011).

An alternative encoding scheme for interval timing would be to directly use the information provided by the tonically active ORNs. Because the information is contained in the temporal profile of these cells, a memory mechanism is needed to transform it to an instantaneously accessible form. One popular method is to implement it through network dynamics; a recurrent network of a set of reliable neurons with short time constant forms a high-dimensional response with a wide range of time constants (Maass et al., 2002; Ganguli et al., 2008; Goldman, 2009; Laje and Buonomano, 2013). However, such implementations are complex and require more biological resources, either in terms of the number of neurons and synapses (Maass et al., 2002), or maintaining precise fine tuning of high SNR (signal-to-noise ratio) neurons (Ganguli et al., 2008; Laje and Buonomano, 2013). In contrast, our proposed encoding scheme is highly efficient, because it requires no connectivity, and uses a small number of noisy sensory neurons.

One of the most robustly observed properties of interval timing across species is the scalar property (Matell and Meck, 2004; Buhusi and Meck, 2005; Jazayeri, 2008). If interval timing is internally kept by counting the cycles of a noisy internal oscillator, the variance of the timing error linearly increase as the target time interval gets longer. However, various experimental paradigms show that the SD linearly increases (known as the scalar property) which means that the timing procedure is not likely to be based on counting. The scalar property has been verified not only in behavior but also in fMRI studies (Buhusi and Meck, 2005). Interestingly, the scalar property can emerge naturally from the bORN population code. We observed that the SD of the interburst interval distribution has a linear relation to the mean (Fig. 5*a*). This gives the main likelihood peak a width that is proportional to the mean, at least to a first order approximation. Therefore, we suggest that such a population of uncoupled oscillators as represented by bORNs could be a source of the scalar property. This is observable from Figure 7. However, a detailed analysis of this possibility awaits future study.

For any sensory input to be meaningful, it needs to be interpreted in a way that ultimately produces meaningful behavior. Although understanding the downstream target(s) and behavioral consequence(s) of bORN input is well beyond the scope of the present study, we can hypothesize that the decoding principle of the bORN population is combinatorial in that, if there were no noise, the population phase distribution of the neurons would be unique up to the least common multiple of the periods, as first noted by Miall (1989). A strikingly similar principle is also observed in the grid cell population in rats (Hafting et al., 2005; Fiete et al., 2008; Sreenivasan and Fiete, 2011), where each grid cell responds in a periodic manner to the relative position of the animal. Then, the downstream neurons (the place cells in hippocampus) decode this combinatorial code and represent the position as a spatial map (for review, see Giocomo et al., 2011). Although the assumption of near perfect periodicity of the code does not hold for the population of bORNs, we suggest that the noise (irregularity of IBIs) inherent to each bORN can be considerably attenuated by encoding the information over the larger population of bORNs, and could be retrieved by an appropriate readout mechanism.

Although the exact downstream neuronal targets (and hence the readout mechanism) of bORNs are unknown, the maximum likelihood (Eq. 8), or maximum a posteriori given a prior distribution of interval distribution of interest, decoding could be implemented within a single projection layer with lateral inhibition implementing a winner-take-all strategy. This possibility is consistent with our understanding of the neural organization of the first synaptic relay in the lobster olfactory brain, the olfactory lobe (OL), thereby providing a neural substrate where decoding potentially could occur. The ORNs project to the OL where they synapse with projection neurons that have extensive lateral inhibitory interactions mediated by local interneurons (Wachowiak et al., 2002). If so, the resulting representation in the decoding layer presumably would be sparse, because only a small population of active neurons represents each time interval. It is biologically plausible, however, that the animal does not directly represent the time interval but rather computes a function of the interval that is most related to a reward signal, with projection weights learned via reward-based synaptic plasticity. Understanding of the readout mechanism, however, awaits future experimentation. Also not addressed by the study is the question of the potential odor specificity of bORNs. Because we used a complex odor mixture to maximize the incidence of bORN responsiveness, we have yet to determine whether bORNs are broadly or narrowly tuned. To the extent that they are narrowly tuned, subpopulations of bORNs could have the added potential of separately encoding the temporal structure of particular odorants. However, one might argue that quality coding could be better done by the much larger population of (more traditionally studied) tonically active ORNs, and bORNs would be more broadly tuned to maximize the probability of stimulus detection for spatiotemporal analysis of the odor plume itself.

In this paper, we have not explicitly modeled possible concentration dependency. We showed earlier (Fig. 4*b*; Bobkov and Ache, 2007) that changes in concentration translates *P _{e}*(ϕ). In our bORN model, this would induce a translation in

*q*(Δ) from Eq. 6 to

*q*(Δ + δ), where δ is the shift in

*P*(ϕ). That means the naive decoder without knowledge of the concentration would shorten (for higher concentrations) or lengthen (for lower concentrations) the estimated intermittency. However, this is not necessarily an inherent confound because presumably local concentration can be estimated from other input, for example, via tonically active ORNs.

_{e}The ensemble coding strategy we report here appears to be specifically adapted to longer stimulus intervals. For most animals, including lobsters, the physiologically meaningful odor signal intermittency ranges from fractions of second to many tens of seconds (Murlis et al., 1992; Vickers et al., 2001; Webster and Weissburg, 2001; Gardiner and Atema, 2010; Reidenbach and Koehl, 2011). It could well be that faster time scales, i.e., fractional seconds, are processed by a different strategy than the bORN ensemble coding strategy we report here, because such short time scale signals match or approach the time scale of typical neuronal dynamics.

Most animals, including lobsters, actively sample their odor environment through a reflex known as sniffing (called flicking in lobsters; Wachowiak, 2011) that in the case of lobsters intermittently moves the olfactory organ through the environment (Koehl et al., 2001) and in the case of vertebrates intermittently flushes the environment past the olfactory organ. Flicking enhances the detection of an odor pulse (Schmitt and Ache, 1979) such that it might be expected to work in concert with the bORN-mediated intermittency. Indeed, this prediction is supported by a recent study of odor plume structure (Reidenbach and Koehl, 2011) showing that the intermittent structure of the odor landscape can be enhanced by flicking, making the relation between the odor encounter interval and the distance to the source relation stronger. Given this sharpening effect, one can envision sampling intermittency (flicking) and the detection of stimulus intermittency evolved to complement each other and strengthen olfactory scene analysis.

The present study represents only a first step toward understanding the roles of rhythmically active uncoupled sensory neurons in information processing. It will be instructive, for example, to compare the activity of bORNs to that of tonically active ORNs, and the response of bORNs to more realistic stimulus intensity profiles to gain deeper insight into the detailed functionality of bORN ensembles. It is also unclear whether different animal species evolved a similar strategy to deal with the temporal characteristics of olfactory signals, such that our findings would have general applicability. As mentioned in the Introduction, both intrinsically and conditional rhythmically active ORNs occur in amphibians and mammals (Sicard, 1986; Frings and Lindemann, 1988; Reisert and Matthews, 2001a,b; Ukhanov, personal communication), and mammalian vomeronasal receptor cells are rhythmically active (Holy et al., 2000; Arnson and Holy, 2011), suggesting the involvement of rhythmically active uncoupled sensory neurons in olfactory processing is not unique to the lobster model and is likely to be a fundamental aspect of olfaction.

## Footnotes

This research was supported by the National Institute on Deafness and Other Communication Disorders through awards DC001655 and DC011859 to B.W.A. We thank Dr K. Ukhanov, Dr D. Hwang, In Jun Park, Y.K. Feng, and A.M. Hein for helpful discussion, comments, and general support.

The authors declare no competing financial interests.

- Correspondence should be addressed to José C. Príncipe, Departments of Electrical and Computer Engineering, and Biomedical Engineering, eb451 Engineering Building, University of Florida, Gainesville, FL 32611. principe{at}cnel.ufl.edu