Spiking Sensory Neurons for Analyzing Electrophysiological Data

Low power consuming biomimetic neurons are considered for use in analyzing electrophysiological data. Starting with a circuit model of a Morris-Lecar inspired spiking neuron, we ﬁ rst investigate the dynamic properties. We demonstrate some of its neuro-computational features including type I and type II excitability, tonic and phasic spiking, spike latency and integration. Electroencephalogram (EEG) signals are then used as excitatory input currents and it is shown that the spiking neurons can provide new insights into brain function. The spike rates of the neurons are employed in a classi ﬁ cation task and shown to yield similar performance compared to one using the frequency dependence. We discuss how this circuit has the potential to signi ﬁ cantly reduce EEG data, improve privacy and lower power consumption for portable EEG systems. © 2020 The Author(s). Published on behalf of The Electrochemical Society by IOP Publishing Limited. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/ by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI

The intersection of neuroscience and computing systems has led to investigations of low power consuming spiking neural networks for use in classification tasks. 1 A wide variety of neuromorphic spiking neuron circuits have been explored, 2 but many implementations are based on the leaky integrate and fire (LIF) model, 3 where a spike train is generated when the input exceeds a threshold. This model is a relatively simple resistor-capacitor circuit and is not able to account for many of the spiking behaviors of biophysical neurons such as spike latency and class 2 type excitability. Other models can capture all of the neuro-computational features but result in more complex electronic circuits. 2 Here, we investigate how a spiking neuron circuit can be exploited for analyzing electrophysiological data, whose origin is the firing of biological neurons. The motivation is to show how the neuroinspired dynamics of the circuit might provide new insights into understanding biological signals and how integrating them into an EEG system has the potential to significantly reduce the energy consumption, reduce data, improve privacy and allow the integration of a classification task on a chip.
We focus on the observation of brain activity through the recording of potentials using extra-cranial electrodes via electroencephalography (EEG) signals. EEG is a non-invasive technique used to understand and monitor brain function, 4 for instance in neurodegenerative diseases, which are increasing with aging populations, investigating and diagnosing epilepsy, 5,6 cognitive psychology and of course brain-machine interfacing. 7,8 The EEG market is growing annually at around 10% per year and one important roadblock is the lack of experts to interpret the EEG results. Thus developing a readout circuit with the possibility of very low power consumption and low overhead that can simplify or directly classify the data is an important objective. The circuit could be integrated into each electrode of an EEG system and directly output significant features of the data, allowing a lower sample rate and possibly a lower bit digital signal compared to conventional systems. It will be able to provide enough information to monitor/classify a given task on a portable device or on chip and thus offers promising perspectives for use in high density neural recording devices. 9 We consider the response of the spiking neuron circuit shown in Fig. 1 that is able to reproduce some biophysical behaviors of real neurons while maintaining an extremely low power consumption and relatively light overhead. This circuit was first introduced in Ref. 10, and was inspired by the dynamical properties of the biophysical model of the Morris-Lecar (M-L) neuron. 11 The M-L neuron describes two channels: a fast ion current (Na + in our circuit) and a slow K + current, and is governed by two differential equations describing the membrane potential V m and the activation variable V k of the K + channel. After describing the circuit, we show how varying the input current or varying the geometrical dimensions of some of the transistors allows the circuit to exhibit different excitabilities, spiking latency, and integration. We then consider how these different operating modes enable the circuit to be used as a "sensor" operating with information encoded by the spike rate of the circuit. We discuss its use in the context of detecting different signatures in EEG data and how this can lead to new insights into their interpretation. Finally, we consider how it can be used in an EEG system and demonstrate how its feature selection capabilities permit classifications comparable with the original EEG data.

The Neuron Circuit
We use the equations governing the analog circuit description found in Ref. 10. Transistors labeled P Na and N K govern the channel dynamics of analogous "Na" and "K" ion channels, with each one connected to a feedback loop. The resting potential of the circuit is at V m < 0 V, preferably near the ∼−80 mV resting potential in the cell. At this voltage, P Na is turned off because the inverter controlling its operation is on and N K , controlled by an amplifier with a gain of 1, is also turned off. When a positive I exc is applied, V m increases and P Na turns on quickly, according to time constant R NA C M , compared to the response at N K , determined by transistors P 2 and N 2 and C K . When C K discharges N K turns on, and causes V m to decrease. This results in a positive inverter output in the controlling circuit of P NA , which will turn it off. With a constant I exc applied, V m will then start to increase again and the cycle is repeated, creating an oscillating voltage (or spike train).
Transistors P Na and N K are thus seen to be analogous to the gating variables in the M-L model, which represent the probability for the Na and K channels to be open. Here they are replaced by the exponential dependence of the transistors the deep subthreshold ds 0 gs t ( ) for an n-type transistor, where G 0 is a pre-factor that includes the transistor's width and length, V gs is the gate to source voltage, V t is the thermal voltage and η is the subthreshold slope factor (ideally h = + C C 1 where C d and C ox are respectively the depletion and oxide capacitances of the transistor). Further details about the correspondence between this circuit and the M-L model can be found in Ref. 12. All transistors in Fig. 1 are assumed to operate in deep subthreshold and using standard circuit theory the equations describing the behavior of the circuit are obtained. As in the M-L model, this results in two differential equations that govern the state variables V m and V k : where: I exc is the excitation current, G Na , G K , G L , G nx , G px correspond to geometrical and device characteristics of the transistors in the circuit, V d , V L are respectively the supply voltage and the gate voltage of the leak transistor and V et = ƞV th . Note that our circuit differs from 10 by the inclusion of a leak transistor, which we will see allows us to tune the threshold of I exc . We model these equations and determine their dynamical parameters using python simulations and the modules Brian2 12 and SymPy. 13 The simulations were calibrated to the circuits realized in a TSMC 65 nm CMOS process that was presented in Ref. 12. To calibrate the simulations, we used C m = 50 fF, C k = 100 fF and I exc = 1 pA to reproduce the simulations in Ref. 12, Fig. 7a, with characteristics in Ref. 12 , Table II. This gives a subthreshold swing S = 92 mV/decade. The calibration parameters result in a very fast spiking neuron and for EEG applications, it is more power efficient to use spike rates that are closer to the observed data. To achieve this, we increased the values of C m and C k , whose absolute magnitude dictates the spike rate of the circuit, via the associated RC time constants. The circuits fabricated in Ref. 12 demonstrated a power consumption of ∼ 60 pW at rest and up to 90 pW and a spike frequency of ∼1.3 kHz. In our circuits we focus on smaller spike frequencies and can expect an even lower power consumption. Compared to the other sources of power consumption in a complete EEG system, as discussed below, this consumption is negligible.
The dynamics of the biophysical M-L neuron have been well explored. 14, 15 The system can exhibit type I excitability where spike trains are the result of a saddle node bifurcation with a spike rate proportional to the magnitude of I exc . It is also possible to observe type II neuron excitability resulting from a Hopf bifurcation where the spike train starts at a finite frequency. In addition, bursting can also be observed, which indicates the presence of multiple equilibrium. The circuit explored here is not a direct one-to-one correspondence with the M-L neuron and does not exhibit all of its characteristics, in particular those involving multiple equilibria. Two coupled neuron circuits, however, could reproduce these properties. While there are electronic implementations that can provide a more exact approximation of the M-L neuron, 16,17 notably by reproducing the shape of the V-nullclines and allowing for several equilibria, the advantage of the circuit in Fig. 1 is the reduced power consumption, the possibility of an integrated circuit realization and the small number of transistors. These properties make it valuable for realworld applications.
To understand the dynamics of the system, the time variation of V k and V m at constant I exc is explored in Fig. 2 with three different sets of parameters, which are reported in Table I. Included in the figures are the nullclines, the curve along which either dV k /dt or dV m /dt is zero, which control the dynamics of the system. The    has only a single equilibria point. As I exc is increased the intersection of the two nullclines occurs at larger values of V k , the limit cycle reduces in size and eventually does not dip below V m = 0, which is our threshold for detecting spikes. The nature of the bifurcation was explored in more detail by computing the eigenvalues of the Jacobian matrix at the equilibria for the different parameters to confirm the nature of the bifurcation (saddle point vs Hopf). Figure 2a demonstrates a type I neuron with tonic spiking starting near I exc = 0 pA, a saddle point potential and the spike rate monotonically increasing with current (Fig. 2b). The parameters in 2a and 2b include a very small G L , which allows spiking at small I exc . A larger G L , as in Figs. 2c and 2d, draws current away from the K-channel and results in a threshold for spiking, even though the neuron remains a type I neuron with a saddle point bifurcation.
A type 2 neuron can be realized by decreasing the gating of the Na channel (G na ) and decreasing the activation of the K channel through G p2 . Changing G p2 through this transition is directly analogous to the one parameter evolution of the M-L neuron that has been investigated in the biophysical model. 15 Figures 2e and 2f show a type II neuron that exhibits damping (sub-critical Hopf bifurcation) before and after tonic spiking. As I exc is increased the neuron exhibits tonic spiking and then becomes sub-critical with the oscillations increasing above V m = 0. Note that at high I exc this neuron exhibits a "phasic-like" spiking, where a single spike is generated and the subsequent damped oscillations are very small. Figure 3 shows the spike rate vs frequency for the neurons from Fig. 2. We observe that the "neuron 1" type excitability exhibits a monotonic increase. For the simulation that used a very small G L , spiking begins near I exc = 0 A and increases from a very small spike rate. For the larger G L , tonic spiking begins later and starts abruptly at a finite spike rate, much like in a type-II excitability. However, no damping is observed and the spike rate dependence on I exc resembles that with a smaller G L but at smaller spike rates. Neuron design for applications would minimize the spike rate in order to reduce power consumption of the device, which implies larger capacitances for C m and C k . Maximizing G L therefore provides an additional parameter that can reduce the spike rate with a comparatively small overhead compared to increasing C m and C k and therefore help reduce the size of the overall neuron circuit.
The spike rate dependence of neuron 2 is quite different from neuron 1. This curve exhibits four regimes: (1) no spiking, (2) damped spiking (sub-critical bifurcation), (3) tonic spiking and (4) damped spiking (super-critical bifurcation), where the spike rate decreases but the damping remains above V m = 0 V. This neuron can therefore provide a much wider range of behavior compared to neuron 1. First, if the input (I exc ) is non-stationary, and the neuron spike rate is approximately the same as the rate of input oscillations, the difference between tonic spiking and damping will be negligible. Essentially, as I exc goes above the threshold for spiking, only a few spikes would be observed before it descends below the spiking  threshold and whether the subsequent spikes are damped or tonic is irrelevant. In the alternative case, the value of I exc reaches the damped region and remains there for many spikes. In this case, the parameters of the neuron can be adjusted so that a decrease in spike rate with increasing I exc is observed, characteristics not possible for a neuron with type 1 excitability. This type of effect, a decrease in spike rate with increasing I exc could be used to identify large artifacts in EEG data. Finally, at very high I exc , the phasic-type response of the neuron can be used as a signal detection-for instance as a detector of an Event Related Potential (ERP) in a brain-machine interface setting. 18 We have now explored the dynamic properties of the circuit with a constant I exc , but EEG data exhibits non-stationary signals and therefore an analysis of the dependence with time is also important. If the time constant of the channels is not too fast, we can observe two related behaviors in this circuit: latent spiking, where a spike occurs after a short pulse, and integration, where a series of pulses can incite a spike, as illustrated in Fig. 4. These effects are generally observable for larger values of C k and C m and for I exc near the beginning of tonic spiking. These types of behavior are the most useful for coupling neurons together, for instance to create frequency filters above a threshold.

Application to Electrophysiological Data
The goal of this section is show how some of the different operating regimes mentioned above can be used to detect different signature of EEG data and thus provide a feature selection of the raw sensor data. We focus on a data set that explores the sleeping patterns of individuals, 19 an important part diagnosing sleep disorders. Recent interest in improving sleep has led to technologies that enable enhanced stage 3 or deep sleep by sending sounds pulses at the properly identified moment. 20 To classify different sleep stages, sleep specialists typically carry out a quantitative analysis of different biophysical tests and different sleep stages are assigned manually to each epoch. 21 We focus on the three stages of Non-Rapid Eye Movement (NREM) sleep: stage 1, transition sleep, stage 2, light sleep, stage 3, deep sleep. Medical practitioners look for particular features for each stage. For instance, stage 3 is characterized by large amplitude EEG patterns compared to the other stage, Figure 5. Sample EEG data for each stage of sleep and for each neuron. The EEG data is in red (left axis), the output from the neurons (right axis) is in blue for neuron 1 and green for neuron 2. stage 2 and 3 exhibit two types of signatures known as K-complex and sleep spindles that distinguish it from stage 1. Recently data scientists have considered this classification problem in the context of machine learning and often use the fast Fourier transform and look for key frequency signatures. 22 The data set we consider was sampled at 128 samples/sec and typically ranges in ± 100 μV, as shown as the red curves in Fig. 5. It consisted of four bipolar EEG signals (C3-M2, O2-M1, C4-M1, O1-M2) and a manually annotated hypnogram was done by a sleep expert, 21 following the AASM rules. 21 The single channel data in Figs. 5-8 used the C3-M2 signal. To simulate the response of the neuron circuit to the EEG data, we assume it is input into the neuron as I exc through a 1 GΩ resistor (normalization of 1 × 10 −12 pA). This normalization can eventually be controlled by a current source using an additional transistor. Figure 5 shows the neuron output, V m , of the two types of neurons (both with G L = 15 pS) for a single 30 s epoch of stages 1, 2 and 3. Clearly visible here is that the spiking appears to increase in stage 2 relative to stage 1 and stage 3 relative to stage 2. This difference is mainly due to the larger amplitude of the signal in each successive stage. We note that the circuit only responds to the positive part of the EEG signal. One can, however, obtain a response for both positive and negative components of the EEG signal through the use of two different neurons.
We observe that both neurons follow the same general trend in spiking patterns but that neuron 2 spikes much less frequently than neuron 1, as might be expected from the plotted spike rates in the 0 to 100 pA range in Fig. 3. We also observe that when the amplitude is too large, as for instance in stage 3 near 1.5 and 19 s, the amplitude for neuron 2 is attenuated, as expected when the neuron becomes damped at higher frequencies. Such attenuation could be useful for eliminating large variations in EEG due to body motions. We have tested the response of the neuron 1 with G L = 200 pS, but not surprisingly there are very few spikes because the threshold to spike is very large compared to the normalization of the current in the 0-100 pA range.
To understand how the different operating regimes might affect the response of the neurons, we reran the simulations using different normalizations. A sample of the responses is shown in Fig. 6 with  the response of Neuron 1 (G L = 15 pS) to three different normalizations of V m in stage 2. Essentially, we see that changing the normalization changes the number of spikes observed. The spike latency is visible near the large amplitudes where many spikes are observed despite the smaller normalization. We have also explored increasing the normalization in a type 2 Neuron to explore the effect of phasic spiking on the behavior. We find that increasing the amplitude results in many more spikes and that it is difficult to isolate phasic spiking by simply changing the normalization. Our investigations suggest that a more detailed analysis involving sequential neuron inputs could greatly improve the detection of individual signatures in the EEG data. Figure 7 shows a hypnogram, where the spike rate for each 30 s epoch is plotted during a night of sleep. Note that in this figure we have removed all data from REM and awake periods of sleep. We observe that there is an increase in the spike rate as the patient moves from transition sleep (stage 1) towards light sleep (stage 2) and then deep sleep (stage 3). However, transitions from more intense sleep towards lighter sleep (stage 3 → stage 2) are characterized by abrupt transitions in the spike rate. In fact, looking more closely, one observes the presence of two slopes in the spike rate vs epoch number as the patient progresses from stage 2 towards stage 3 sleep. There is a first sharp slope over several minutes as the patient moves towards stage 3 sleep and then a continued increase in the spike rate but with a smaller slope once the patient is in deep sleep. The pattern of these changing slopes evolves over the night of sleep, in particular developing a prominent sawtooth pattern as the night progresses. This unusual way of exploring EEG data suggests that it could provide new understanding of electrophysiological data.
Finally, we consider using the spike rates in the classification of the 3 stages of sleep. We consider that the spike rate determined for each 30 s epoch of EEG data, will provide sufficient information to detect a given stage of sleep. We compare the output of a single channel of EEG for both neuron 1 and neuron 2 separately and then together. We also used the output of neuron 1 for 4 different EEG electrodes for an additional classification. Using the sci-kit learn python library, 23 we employed three standard machine learning classifiers: random forest, support vector machine and Naïve Bayes. The first two have been widely explored for automatic sleep detection. 22 The Naïve Bayes classifier is highly sensitive to the independence of the features and has been shown to give very good performance for sleep stage detection when few feature are used. 24 As a comparison, we also classify the results by taking the Fast Fourier Transform (FFT) of the EEG signals and using as inputs the power spectral density of the frequencies between 0-20 Hz. Because the stages are unequal in size, the training data is padded so that each stage has the same number of points. The test results however are classified using the same proportion as the initial data. We used 70% of the data set for training and 30% for testing. We performed each classification 550 times with different random selections of the data to obtain an average accuracy. The average results and standard deviations are summarized in Table II. The main result is that the spike rate alone can give as good a classification as the FFT of the EEG data, although not necessarily with the same classification technique. In particular, Naïve Bayes gives the best classification for the spike rates. This is not surprising given that it requires the inputs to be linearly independent and we have used very few features (just a single point for the classification of each epoch for a single EEG Figure 8. Comparison of the confusion matrices for the best classifiers of the FFT data (random forest) and the spike rates (Naïve Bayes) for a single EEG channel. Data from the same training and test sets were used. The data in the single channel classification was from the EEG C3-M2 signal. For the 4-channel classification, for the FFT data, an average of the 4 signals (C3-M2, O2-M1, C4-M1, O1-M2) was used and for the spike rates the individual outputs of each neuron was used as a feature. The FFT data was obtained for each 3840 datapoints (30 s) of the raw EEG data. The classifications involving spike rates involved just a single rate at each 30 s. channel). We also see that the additional EEG channels increases the accuracy of the FFT classification more than that of the spike rates classification. To compare the individual categories of the best classification using the FFT and the spike rates, in Fig. 8 we show the confusion matrix for a Random Forest classification of the FFT data and the Naïve Bayes classification of the spike rates. We see that the spike rate classification using Naïve Bayes does best at classifying stage 3, mainly because the large amplitudes result in a higher spike rate and worst distinguishing stage 1 and stage 2. The FFT classification using the Random Forest is best at classifying stage 2 and then in general performs better on stage 1 than the spike rate classification.

Discussion
The sensor readout circuit proposed here is useful for wearable EEG systems 25,26 or high density neural probes. 27 Tt acts as a feature selection and can significantly reduce the data and therefore the energy consumption. The main sources of power consumption in portable EEG systems (that do not have on-chip classification) are 25 : (a) the amplifier, (b) the analog to digital converter and (c) the data transmission to a device in close proximity such as a cell phone or computer, that can either perform a classification and/or send the data to the could for further analysis. A typical state of the art low power consumption amplifier uses ∼2 μW per channel. 28 Commercially available analog to digital converters such as the T.I. ADS129x chips, which have multi-channel delta-sigma converters as well as programmable gain amplifiers, typically consume 0.75 mW/channel. Recently published research has demonstrated a 15 bit delta-sigma converters with a 20 μW power consumption. 29 The main bottleneck in portable EEG systems is the third component: energy consumption of the data transmission, which is proportional to Jf s R where J is the energy per bit (∼ 50 nJ/bit), f s is the sample rate and R is the resolution, typically 12 bits. To surmount this, researchers have begun to move towards onchip classification. Recent research explored on chip classification with 9 RISC-V cores and a microcontroller unit to demonstrate a power consumption of 6.31 mW for an 8 channel system. 30 Another recent paper showed an on-chip decision classification with a 71 mW power consumption for 4 channels but just 0.9 mW for the readout electronics.
The readout circuit of EEG sensor data has two potential uses. The first is to consider the spiking neurons as a feature selection that could significantly reduce the energy consumption during data transmission. In this case, the readout circuit still requires the amplifier but the spiking neuron acts as the ADC. It has been shown to have ∼40 fJ/ spike of dynamic power consumption 10 and if we consider the results of the spike rates in Fig. 7, this would imply ∼ < 100 pW of power per electrode. It does, however, require the spikes to be converted into a rate, which would dominate the power consumption. This could be done using either a circuit involving a floating gate transistor, 31 where the power consumption can be as low as 10 nW, 32 or one of the many novel emerging memory devices, 33 which would also allow the power consumption to remain very low. The data resolution of the spike rate is limited by the number of possible spikes in the time period and for neuron 1 we observed that a 7 bit resolution would suffice. More importantly, the sampling rate of the data transmission is substantially reduced because only a single value per channel and epoch would be transmitted compared to the 3840 values of raw EEG data obtained for a sampling rate of 128 samples/sec. The bottle neck of energy consumption would therefore no longer be dominated by the data transmission (on order of 100 s of μW) but by the ∼ 2 μW amplifier. Note that while we have used data from traditional EEG sensors here, because the input of our readout circuit involves an excitation current, a promising direction to reduce the power consumption of the amplifier might be to use the transducer realized with organic electrochemical transistors, 34,35 which output a current directly into the spiking circuit.
The circuit can also pave the way towards a very low energy classification on-chip. A very simple implementation for instance could consider the spiking neurons as stochastic bit-streams and employ the C-Muller elements to implement a Bayesian inference. 24,36 Such an implementation would first require a supervised learning step for different patients so that the priors for the different classes of data could be obtained and subsequently implemented in the circuit using a spiking neuron at fixed I exc . An alternative method for a more complex sensor network involving many features would be to implement the testing on chip in the form of a spiking neural network but to realize the training offline by converting the spiking network into an artificial neural network that can be trained offline. 37 More generally, this type of circuit could potentially be used to record any biological non-stationary electrical signal such as intracranial electrodes (ECoG), muscles using electromyography (EMG), 38 and the heart via electrocardiography (ECG). 39 To realize such different applications, the parameters of the neuron circuit can be engineered to respond to the important frequencies and signatures of the different measurements. The broad applicability of this technique to a large array of different biological sensors suggests that using spiking neuron circuits could allow facile sensor fusion in an integrated system that used a spiking neuromorphic architecture for classification of a given task.

Conclusions
We have explored the dynamics of a spiking sub-threshold neuron circuit. The circuit is able to reproduce the dynamics of the two types of neuron excitabilities. We explored using the circuit to probe the non-stationary time series signals of electroencephalography. Our results show that the average spiking rate of the neurons can capture the amplitude characteristics of this data, resulting in a classification precision comparable to spectral classifications using conventional machine learning techniques. We also investigated the impact of different neuron characteristics. Such manifestations could be used to realize parts of the EEG detection, for instance reduction of large amplitude artifacts, analog to digital conversion or distinguishing large amplitude events. The manifestation of electrophysiological signals is due to the temporal synchronization of neural activity. Using sensory circuits inspired by modeling the behavior of neural activity to probe such signals thus provides natural synergies that should enable enhanced detection with much lower power consumption. Such system integration holds a potential for groundbreaking innovations in ultra-low-power biomedical circuit design and flexible bio-sensing circuits.