Classifying continuous, real-time e-nose sensor data using a bio-inspired spiking network modelled on the insect olfactory system

Chemosensing electronic nose ('e-nose') technology has great potential for applications in everyday life, ranging from drug detection to food quality assessment [1] and even diagnosis of illness [2, 3]. A typical e-nose device might employ doped tin oxide sensors that generate time series of resistance values that change selectively in the presence of certain chemicals. Typically, pre-processing and encoding methods are applied to convert the collected time series into discrete numerical representations of each odour that are amenable to separation by classical classifier algorithms such as k-means clustering or support vector machines (SVMs) [4]. However, e-nose devices have had limited success in some applications because odour detection and classification comprise a highly challenging domain characterised by high dimensionality, unknown organisation of the vast 'odour space' of all volatile chemicals and complex dynamics of odour plumes [5]. To compound these difficulties, current sensor technology continues to exhibit distinct shortcomings in speed, sensitivity, selectivity, recovery, and drift avoidance [3]. As a result, e-nose devices remain substantially outperformed in both speed and accuracy by biological noses, both mammalian and insect. This is particularly true in realistic field settings outside the controlled environment of the laboratory. Whilst development of novel bio-based sensors has begun to show promise [6–8] we focus here on another challenge, the classification of continuous incoming olfactory data, leveraging the intrinsic temporal nature of spiking network approaches. In the work presented here, we use a benchmark data set of 20 chemicals measured by 12 metal oxide e-nose sensors [4] to investigate an odour classifier employing a spiking network design abstracted from the insect olfactory system.

Olfactory regions in the insect CNS (figure 1(a)) follow a common organisation, comprising three distinct stages of sensing, transformation and association (classifying) [9]. At the sensing stage, neural activity (spiking) is triggered in olfactory receptor neurons (ORNs) located along the antennae. Each ORN expresses one type of olfactory receptor (OR) protein which binds a different and specific range of ligands that may appear in a variety of airborne chemicals [5, 9, 10]. The net effect is that each odour will trigger a fast, widespread spiking response across the set of ORNs, with varying intensity in each ORN type, depending on the odour's chemical constituents. It is likely that this distribution of response profiles has been tuned through evolution to best aid discrimination between odours.

Zoom In Zoom Out Reset image size

The transformation stage takes place in the antennal lobe (AL) structure where axons from the ORNs synapse onto a limited set (e.g. 43 in Drosophila) of spatially distinct sites of high synaptic interaction called glomeruli [11, 12]. ORNs that express the same type of receptor [9] provide excitatory input to subpopulations of projection neurons (PNs) and local interneurons (LNs) within a given glomerulus. PNs generate the output from the AL to higher brain centres. LNs, in contrast, inhibit PNs in other glomeruli in the AL, shaping patterns of glomerular activation as they compete. Whilst the exact coding scheme enacted by this competition remains contentious [13–15] it is evident that the broadly tuned ORN responses are modulated by this lateral inhibition to generate a sparser and more granular encoding intended to enable separation of even closely related odours.

Activity from the AL is projected to higher brain centres like the mushroom body (MB) and the lateral horn, where classification and multisensory integration take place [16]. Here, the activation patterns of the glomeruli are ultimately associated with an output that is relevant (useful) to the insect, for example the presence of a foodstuff.

A number of neuronal network classifiers based on the insect olfactory system have been developed and have demonstrated some success in classifying static datasets comprising either standard benchmark sets [14, 17] or artificial data [18]. However, these models remain largely unproven in a practical setting where input sensor data is real-time, continuous, potentially noisy, lacks a precise onset signal and accurate classification requires the inclusion of temporal aspects into the feature set. This investigation seeks to inform and develop the potential and suitability of biomimetic classifiers for use with typical real-world e-nose sensor data.

Previous classifying work with this data set [4] has shown that downsampling the data but retaining temporal features in the encoding process produces compact representations of the data that allow successful recognition of odours by a linear SVM classifier. This approach also out-performed a range of non-temporal encodings [4]. For example, creating a single 16-dimensional vector representation by sampling a response at only four time points on just 4 sensors (figure 2(b)) was sufficient to classify the dataset to more than 99% accuracy using the SVM. Furthermore, restricting the sample times to the early stages of response—for example the first 30 s out of 300 s total—caused very little reduction in performance [4]. Based on these results, we chose a similar encoding approach to create a static multivariate representation as input to our own bio-inspired neural model, namely, four sensors, sampled at four time points from odour onset (figure 2(b)).

We generated our set of VR points to represent the 16-dimensional input space by using a self-organising 'neural gas' process [19] against the complete training dataset. Note that this is unsupervised, i.e. makes no use of class labelling. A representative example of the mapping results is provided in figure 3(a).

Zoom In Zoom Out Reset image size

One of the practical problems introduced by this timed-sample encoding strategy is that, in the field, an e-nose device will rarely have access to accurate odour onset timing. The requirement for acquiring precisely timed temporal features from data that is continuously varying without timing cues therefore appears particularly problematic.

However, by adding delay lines to the input, it becomes possible to access time-shifted versions of the input signal which correspond to the required intervals between the chosen sample points (figure 2(c)). When the input signals are read the delayed signals are also accessed from the delay line, generating a compound input vector comprising both current and past samples (figure 2(c)). If the employed delays used are the same, then at some point in time, as the sensors' response unfolds, this compound vector will reproduce the same timed sample set used to train the classifier.

This approach is not, we suggest, a purely engineering convenience that diverges from the bio-inspired basis, as a conceptually similar approach to capturing temporal features has also been proposed to form part of the insect strategy for odour classification [20]. For example, the axonal 'delay-line' effect introduced by the spaced layout of the receptor array in the Acheta Domesticus cricket has been shown to modulate the neural coding of odours [20]. Although these delays are clearly far shorter than the sample intervals employed here, it should be noted that the sensor response of insects are orders of magnitude faster than our metal oxide sensors and would indeed require very fast sampling to implement temporal feature encoding.

Although this delay-line strategy can address the sample timing issue with continuous data, there is also the problem of knowing when to appropriately begin and end evaluation of spiking activity at the output. This scenario can clearly occur in a detector or alarm deployment where sensors are left open indefinitely, awaiting the presence of a target odour. To address this, we tested a simple spiking activity-based algorithm, which identifies sufficiently active periods and demarcates appropriate evaluation time windows over which a comparative spike count can be made (figure 4). This allows the classifier to function without any use or knowledge of a global or recording-based time variable.

Zoom In Zoom Out Reset image size

The classifier was trained and tested using a stratified ten-fold cross validation strategy. On each of ten runs, one (different) member of every class is allocated to the test set. The remaining 180 observations (recordings) are randomly drawn for presentation from each class in turn, with random class ordering per run.

The response to the temporal signature (static representation) of the selected recording is applied as input to the Poisson 'receptor' neurons and, with plasticity enabled, the spiking network simulation is run for a (simulated) 500 ms using an integration timestep of 0.5 ms.

Exposure to the training set is repeated a fixed N times until performance on the training set has converged, without straying into overfitting (performance on test set drops). The value of N was set during initial investigations—see figure 6(d).

Without using a correlate of the widefield receptor response, the AL model acting directly on a stratified 20 class dataset (encoded as RN spiking activity driven by the normalised response per sensor) was able to classify a test dataset to an average 39% (s.d. 2.7%) across ten-fold cross-validation. This suggests that the classifier is functional but that the dataset comprises a non-trivial classification problem that the relatively simple strategy enacted by the AL structure alone does not solve easily.

The input was switched to the proximity response to a set of VR points in input space generated independently from each training dataset in the cross validation. All further results were elucidated employing this two stage model.

Figures 5(a) and (b) provides a representative example of the spiking activity and cluster interaction of the classifier when presented for 500 ms with the static encoding of an unseen sample of isopentylacetate (odour class 13), after training is completed. The pattern of spiking activity in the RN clusters demonstrates how an input elicits a wide ranging response encompassing a large proportion of the VR receptors. The 2D projection of the response (figure 5(c)) illustrates how the response is proportional to the proximity to the input. The activity in the PN clusters demonstrates the action of the inter-cluster mutual inhibition to reduce the activity to the closest VRs whilst emphasising the difference between these. Finally, the AN cluster activity illustrates how the correct class cluster has been associated, via the appropriate strengthening of PN–AN connections during training, with the VRs responding most to samples from this class. The lack of activity in rival AN class clusters also reflects the effect of the WTA inhibition in this layer.

Zoom In Zoom Out Reset image size

After training, the classifier performance was tested in three scenarios of increasing difficulty. Firstly, plasticity is disabled and, as with the training, the static representations of the test set recordings are simply applied to the classifier in turn for 500 ms. The classifier's winner decision for each recording is awarded to the AN class cluster with the maximum spike count over the 500 ms presentation. This simpler scenario was used primarily for parameter tuning and to characterise the classifier performance.

In the second scenario, performance is tested by replaying the continuous timeseries recording data from the e-nose. The sensors' response was originally captured at 2 Hz (i.e. 500 ms between samples) totalling 600 samples over a 300 s recording. The recording is played into a software delay line to generate three additional delayed versions at 7 s intervals. From these four inputs sampled together the VR response is obtained and applied to the spiking network input for 500 ms. The winner is awarded to the most active class cluster over the whole 300 s presentation.

Finally, performance is tested by adhering to a more realistic 'detector' scenario where there are no preset discrete periods to evaluate over and a positive odour detection (i.e. classifying decision taken) is decided solely on the basis of sufficient activity in the network taking place.

To achieve this, the current average spiking rate of each cluster is calculated on-the-fly and an algorithm is applied to trigger the start and end of an appropriate 'evaluation window' in the time domain (figure 4). The winner is evaluated from the most active class cluster within this window.

Figure 6 characterises the performance of the classifier in the initial testing scenario. Three important configuration settings are explored which impact the resource and processing load required by the classifier. Firstly, the number of VRs employed (figure 6(a)) affects the resolution of the input space mapping but each one adds two extra neuron clusters and associated synapses to the model. Results suggest that, for this data set distribution and self-organising algorithm combination, little performance is gained beyond 45 VRs.

Zoom In Zoom Out Reset image size

Secondly, the use of probability-based rate coding of the input means that the cluster size (figure 7(b)) affects the reliability and consistency of the mutual inhibition behaviour that is an important part of the classifier strategy. This comes about through issues such as synchrony artifacts caused by the resolution of the integration timestep (see [14] for a detailed discussion). Results suggest that 30 Poisson neurons is a minimum cluster size and that 50–60 neurons can improve performance further by a few percentage (90 up to 94.55%).

Thirdly, it is clearly advantageous in many situations for an e-nose device to make an identification as soon as possible. Whilst the sensor recordings show that their response does not settle for at least 2 min after onset, results suggest (figure 6(c)) that taking four samples in the period from 10 to 30 s is sufficient for a strong classification performance and relatively little improvement is available by widening this time window.

Finally, performance can be improved by re-exposing the training set a number of times as the association can continue to be refined by the learning rule. However, more repetitions slow the learning process and more importantly, risk over-fitting to the training data over the test data. Results suggest (figure 6(d)) that little performance is gained beyond three repetitions, and that over-fitting sets in at six repeats.

Overall, we conclude that, drawing upon only the first 30 s of this 200 observation, 20-class data set, the described classifier can perform creditably whilst still lagging a SVM (we achieved 94.55% avg. accuracy, 10 × 10-fold stratified cross-validation, against 99% for SVM) when using 43 VRs, a cluster size of 60 neurons and four exposures to the training data.

Considering the question of why our classifier is outperformed by an SVM approach, we examine the specifics of the error cases. We find that 50% of errors are accounted for by the three most common class pairs that are confused when using 30, 43, or 50 VRs (Z2-hexen-1-ol and furfaral, acetone and 2-heptanone, pentanol and hexanal). These can all be shown to involve overlapping or closely positioned clusters in input space (see the figure 4(a) PCA plot for a visualisation).

Figure 7 provides a representative example illustrating the spiking activity in the network when presented with the first 40 s of a 300 s previously unseen continuous timeseries recording of the sensors' responses. As discussed in the methods, the network response is elicited after training has been completed using 500 ms presentations of the static downsampled encodings of the training set, using four samples of four sensors taken at 9, 16, 23 and 30 s after odour onset, (i.e. with a 7 s inter-sample interval).

Zoom In Zoom Out Reset image size

Following training, an unseen continuous recording is tested as input, using a delay line to recreate the sample intervals used in training, and the 4 × 4 result is played as input into the classifier network. From figure 7 it can be clearly observed that the response of the network builds to a peak at the 30 s mark before dying away. This peak corresponds, as expected, to the point where the four input copies should align to the last of the four sample times (30 s) used as the input space for the distribution of VRs and where training of VR-to-class association has taken place. In particular, it can be noted from the heatmap (figure 7(b)) the intensity of response in the AN layer indicating a dominant class cluster. If the 'winning' class is inferred from the AN cluster spike counts over the full available 300 s presentation, this results in a 92% (std.dev. 0.68% across 10 × 10-fold cross-validation assay) accurate classification of the test set.

Disallowing use of these discrete 300 s windows for spike counting make for a more realistic scenario without onset information. In this scenario, windows of high activity from the network are used to infer that a valid classification is being made. The red square in figure 7(b) indicates the evaluation window chosen by the algorithm described in figure 4 and demonstrates the identification of a zone of interest without use of absolute time or onset information. If the 'winning' class is now inferred by the AN cluster spike counts within these evaluation windows, this results in the same 92% classification of the test set (std.dev. 0.56% across 10 × 10-fold cross-validation assay). This result compares favourably with 94.55% achieved using static encoding of the test set, which requires precise onset timing as a pre-requisite.

A final test to validate any classifier and discount effects of inadvertent overfitting due to multiple re-use of data when developing the model is to expose it to an entirely unseen data set and to report the results with no further parameter modifications beyond retraining with the new training data set [21]. Applying the classifier to a separate but directly comparable data set, which was produced one year later using the same procedures resulted in 92.65% (0.74% std.dev., 10 × 10-fold crossvalidation) accuracy for static encoding and 89% (0.61% std.dev., 10 × 10-fold crossvalidation) when continuous input was supplied and the evaluation window was applied. The small but significant drop in the mean performance (confirmed by t-test) suggests that, within bounds of expectation, the classifier generalises well to unseen data.

Data here refers to the GPU-based implementation (figure 1(c)) of the conceptual classifier model (figure 1(b)).

We use clusters of 60 neurons throughout to ensure consistent averaged spike rate coding and to minimise potential synchrony effects occurring in inter-cluster competition [14]. Using 40 VRs, this implies 6000 neurons and some 6 million active synapses. To achieve real time classification at this scale a GPU-accelerated neural simulator (GeNN, [26, 27]) was chosen, using CUDA C to run at 8 times real-time on a standard NVidia card (GeForce GTX 760–1152 cores). Classifier code is implemented in C++ running on a high end Linux workstation (8-core, 3.7 Ghz Xeon, 32GB RAM) with the g++ compiler. Both model and classifier are available as examples in the GeNN repository [26].

E-nose recordings comprise 300 s, 12-sensor timeseries data at 2 Hz sampling. Recordings are zeroed when the chemical is released. A basic smoothing filter is applied using a simple moving average. Timeseries data is then reduced to a static representation comprising a single 16-dimensional vector constructed from four sensors, sampled at four points. Sample time and sensor choice were informed by previous classification work on this data set [4], resulting in a choice of two doped tin oxide sensors and two novel chromium titanium oxide with sample times set at 9, 16, 23, 30 s (see figure 2(b)).

With each recording in the training set now represented by a point in 16-dimensional input space, a set of VR points is selected to representatively span this space, generated for each cross validation dataset using the 'neural gas' algorithm [19].

A proximity response r $\in$ (0,1] is elicited from each VR according to its distance d in data space to the current input point. The response r is made unit free by scaling with the average distance d_avg between all input points in the set. The wide distribution of relatively tight class clusters in data space (figure 4(a)) means that using a linear response function (see [14]) performs poorly on this data set. An exponential decay function is employed instead. The parameters k and m are used to widen the response field whilst emphasising the nearest VR (figure 3(b)). The rate for the ith VR is given by

$$\begin{eqnarray}&&{r}_{i}=\mathrm{exp}\left(-{\left(\displaystyle \frac{k{d}_{i}}{{d}_{{\rm{avg}}}}\right)}^{m}\right),\end{eqnarray} \tag{ 1 }$$

where d is calculated as the 'Manhattan' distance between points, i.e. the sum of the absolute coordinate differences. The response curve is shown in figure 3(b), for k = 5, m = 0.7.

We model an abstracted insect olfactory system in GeNN [26] as three main 'layers': RNs, PNs and ANs as shown in figure 1(c). Sub-population 'clusters' represent VR activation (RN), glomeruli (PN) and recognised classes (AN).

The RN input neurons implement Poisson rate coding, producing spike trains with an average spike rate that tracks the current response level of their assigned VR. The required average spike rate ρ is set as:-

$$\begin{eqnarray}&&\rho =r({\rho }_{{\rm{max}}}-\;{\rho }_{{\rm{min}}})+\;{\rho }_{{\rm{min}}},\end{eqnarray} \tag{ 2 }$$

where r is the (0,1] response of the associated VR and ρ_max = 70 Hz, ρ_min = 5 Hz are the maximum and minimum spike rates.

PN and AN layers are implemented as 'map' neurons [28]. This model employs a phenomenological, discrete time dynamical map with a fixed timestep Δt set at 0.5 ms. These can be computed much faster than conductance models (HH) enabling modelling of up to 10⁵ neurons in real time.

Synapses are considered to have zero axonal transmission delay and the synaptic current I_syn at time t flowing from a presynaptic to a postsynaptic neuron is defined as:-

$$\begin{eqnarray}&&{I}_{{\rm{syn}}}(t)={g}_{{\rm{syn}}}S(t)({V}_{{\rm{rev}}}-{V}_{{\rm{post}}}(t)),\end{eqnarray} \tag{ 3 }$$

where g_syn is the synapse conductance set via the connectivity matrix and is potentially updateable through plasticity. S(t) represents the amount of neurotransmitter available at the postsynaptic neuron. Each presynaptic spike delivers a fixed amount which is decreased exponentially on each timestep. V_post(t) represents the membrane potential of the postsynaptic neuron at time t whilst V_rev denotes the synapse reversal potential (V_rev fixed to 0 mV for excitatory, −92 mV inhibitory).

Topology comprises RN–PN connectivity, PN–PN, PN–AN and AN–AN (see figure 1). Note that all connectivity is inter-cluster, there are no internal connections between neurons of a single cluster. Fed by a VR, each RN cluster excites its corresponding PN cluster (a notional 'glomerulus') using a common fixed weight with 50% random connectivity.

Competition between glomeruli is implemented by PN–PN inhibitory synapses. PN neurons connect only to those in other glomeruli (random 50% connectivity) using a fixed weight chosen to damp activity outside of the most active clusters.

Associating VR activation with output classes is controlled by the PN–AN synapse population, set at 50% random connectivity. Weights are initialised and modified in training according to the plasticity rule (see next section).

The final synapse population, AN–AN acts to provide a single clear 'winning' class cluster in the AN output layer. A WTA configuration is constructed by lateral inhibition between neuron clusters representing each class. In contrast to lateral inhibition between PNs in different glomeruli the inhibitory weight is set higher in the WTA to eliminate activity outside the winning cluster.

Before training runs all PN–AN synapses are initialised randomly between minimum and maximum weights, w_max and w_min. Simple reinforcement style learning is applied to strengthen PN–AN connections that lead to correct classifications and weaken those producing incorrect classifications. To this end, after each recording presentation, weights are adjusted up or down according to spikes in the 'winning' AN class cluster. For the set S_win comprising every AN neuron that spiked in that cluster, all incoming synapse connections from PN are considered and its weight is adjusted by a fixed amount ±Δw, (w bounded to w_{min, max}). The sign of the adjustment follows whether the classification decision was correct or not. Thus the weight w_ij from the ith PN neuron to the jth AN neuron is adjusted as:-

$$\begin{eqnarray}{w}_{ij}=\left\{\begin{array}{cc}{w}_{ij}+{\rm{\Delta }}w, & \begin{array}{cccc}j\in {S}_{{\rm{win}}}, & {\rm{winner}}={\rm{correct}}, & {\rho }_{i}\gt {\rho }_{{\rm{thresh}}}, & {w}_{ij}\lt {w}_{max},\end{array}\\ {w}_{ij}-{\rm{\Delta }}w, & \begin{array}{cccc}j\in {S}_{{\rm{win}}}, & {\rm{winner}}\ne {\rm{correct}}, & {\rho }_{i}\gt {\rho }_{{\rm{thresh}}}, & {w}_{ij}\gt {w}_{min}.\end{array}\end{array}\right.\end{eqnarray} \tag{ 4 }$$