NET-TEN: a silicon neuromorphic network for low-latency detection of seizures in local field potentials

Horizontal hippocampus-cortex (CTX) slices, 400 µm thick, were prepared from male CD1 mice 4-8 weeks old. Epileptiform activity was induced by treatment of the K+ channel blocker 4AP (250 µM). A detailed description of the methods can be found in [38]. All procedures have been approved by the Institutional Animal Welfare Body and by the Italian Ministry of Health (authorization 860/2015-PR), in accordance with the National Legislation (D.Lgs. 26/2014) and the European Directive 2010/63/EU. All efforts were made to minimize the number of animals used and their suffering.

Extracellular local field potentials (LFPs) were acquired using the Mc_Rack software through a $6\times 10$ planar MEA (Ti-iR electrodes, diameter 30 µm, inter-electrode distance 500 µm, impedance ${\lt}100$ kΩ) connected to a MEA1060 amplifier (all from Multichannel Systems, Germany). The brain slices were continuously perfused at ${\sim}1$ ml min⁻¹ with artificial cerebrospinal fluid containing 4AP (see [38] for composition), equilibrated with carbogen, and maintained at 32 ^∘C. The signals were sampled at 2 kHz and low-passed at half the sampling frequency before digitization.

Discrete-time LFP data were encoded in MATLAB using the step-forward encoding (SFE) algorithm defined in [39]. The encoding process gives rise to two spike trains: UP spikes, associated with positive-going signal deflections, and DW spikes, associated with negative-going signal deflections. For each 100 s long LFP sample, two arrays of length 100 s*2 kHz were declared and their elements were initialized to '0'. One array was allocated to store UP spikes and the other for DW spikes. The algorithm was initialized by taking the first value of the LFP signal as baseline. At every new time step, the signal was compared to the baseline ± a threshold. When the value exceeded the upper limit of the interval an UP spike was registered by assigning a '1' to the corresponding element of the corresponding vector and the baseline was updated to the baseline plus the threshold. Likewise, a DW spike was recorded if the lower limit was surpassed and the baseline was shifted to baseline minus threshold. The threshold value determines the sensitivity of the encoder to changes in the signal and the resulting spiking frequency, with smaller values leading to higher spike densities. Here it was selected following an iterative process as the one that returned an output with an average firing rate of approximately $30\%$, where the output was the logic sum of UP and DW sequences. A second script in Python converted the two digital streams into real waveforms. Each '1' was replaced by a pulse of 130 µs full width at half maximum and 200 mV amplitude and the rest was zero padded. Finally, said waveforms were resampled at 10 kHz, saved as two separate columns in a text file and uploaded on a Keysight 33 600 A arbitrary waveform generator.

Silicon neurons were based on the mathematical model of spiking neuron developed by Izhikevich [40]. The model is derived from the bifurcation analysis of the Hodgkin–Huxley model and retains its most salient features, thereby enabling the faithful reproduction of a wide range of firing patterns commonly observed in biological cortical neurons. Izhikevich neurons are described by the following two-dimensional system of ordinary differential equations:

$$\begin{align} \left\{\begin{array}{lr} \dfrac{dV_\mathrm{mem}}{dt}=0.04{V_\mathrm{mem}}^{2}+5V_\mathrm{mem}+140-U+{I}_\mathrm{syn} & \\ & \hfill (1a)\\ \dfrac{dU}{dt}=a(bV_\mathrm{mem}-U) & \hfill(1b) \end{array}\right. \end{align}$$

with the after-spike reset defined as

$$\begin{equation} \text{if}\;\;\; V_\mathrm{mem} \geqslant 30 mV: \begin{cases} & V_\mathrm{mem}\ \text{is set to } c \\ & U\ \text{is incremented by } d \end{cases} \end{equation} \tag{ 2 }$$

where $V_\mathrm{mem}$ is the transmembrane voltage, that is the difference in electric potential between the interior and the exterior of the cell membrane, and U is a state variable named the slow membrane recovery variable. Unlike the Hodgkin–Huxley neuron, the model does not explicitly characterize the electrodiffusion kinetics of single specific ionic channels. However, the slow variable U reflects the overall state of activation of the potassium channels and inactivation of the sodium channels [40]. According to (1a), U acts as a negative feedback for $V_\mathrm{mem}$. In equation (1b), a can be regarded as a time constant which determines the rate of evolution of U, and b is another parameter that defines the strength of the coupling between $V_\mathrm{mem}$ and U, meaning how much the variation of the transmembrane voltage influences the dynamics of the slow variable. Instead, c and d from equation (2) are the after-spike reset values for the two state variables [40].

Synapses adhered to the spike-timing-dependent plasticity (STDP) rule, an unsupervised learning algorithm modeled after biological evidence [41]. According to STDP, synaptic modification depends on the temporal relationship between pre- and post-synaptic firing. This dependency holds true only if the distance between pre- and post-synaptic spike is within an interval of tens of milliseconds. This algorithm provides that, when a pre-synaptic spike anticipates a post-synaptic spike, the synaptic efficacy is strengthened; on the contrary, when the pre-synaptic neuron fires after the post-synaptic neuron, the connection between them is weakened [42]. In mathematical terms, STDP is modeled as:

$$\begin{equation} W(\Delta t) = \begin{cases} A_+ e^{-\frac{\Delta t}{\tau_+}} &\text{if}\ \Delta t \geqslant 0\\ -A_- e^{\frac{\Delta t}{\tau_-}} &\text{if}\ \Delta t < 0 \end{cases} \end{equation} \tag{ 3 }$$

where $W(\Delta t)$ is the synaptic modification; $\Delta t$ is the time elapsed from when the pre-synaptic neuron emits a spike to when the post-synaptic neuron spikes; $\tau_+$ and $\tau_-$ are time constants in the order of 10 ms, and the parameters $A_+$ and $A_-$ are hyper-parameters that set the synaptic weight change rate.

All CMOS circuits were designed and simulated in sub-threshold analog domain by means of Cadence^® Virtuoso^® platform. Full custom layout design of the chip was carried out using Virtuoso^® Layout Suite. Calibre was used to perform the Design Rule Check and verify whether the geometric constraints imposed by the semiconductor manufacturer were met. The same tool was employed for the Layout Versus Schematic verification and to extract the parasitic resistance and capacitance. Parasitics were then harnessed in post-layout simulations to ensure that each instantiated physical cell and its corresponding schematic circuit were functionally equivalent. The chip was fabricated by Taiwan Semiconductor Manufacturing Company in 180 nm process technology. The silicon die was packaged with a 48-pin Quad Flat No-Lead (QFN48) package. The ten output channels of the neuromorphic device were buffered on the printed circuit board (PCB) to increase their drivability. For this purpose, we chose the MAX44280 operational amplifier by virtue of its 0.4 pF input capacitance that limited the load at the output nodes. To minimize the parasitic capacitance introduced by the PCB traces, in the PCB layout, the buffers were placed as close as possible to the neuromorphic chip footprint. Circuit biases were controlled through adjustable voltage dividers on the PCB obtained from 100 kOhm potentiometers with 25 turns to allow fine adjustments. The maximum achievable value for each voltage bias was 200 mV, with one turn corresponding to a variation of 8 mV. $V_\mathrm{DD}$ and $V_\mathrm{D}$ could be tuned up to 450 mV, at the expense of the tuning precision, that dropped to 18 mV per turn.

To find the precision of the system, each measurement of NET-TEN output firing was compared with the label signal associated with the delivered input sample. Precision is defined as:

$$\begin{equation} \mathrm{Precision} = \mathrm{\frac{TP}{TP+FP}} \end{equation} \tag{ 4 }$$

where TP means true positives and it is the number of spikes emitted by NET-TEN in correspondence of any pathological activity, that is when the label was either 1 (interictal) or 2 (ictal); instead, FP are the false positives, the number of spikes registered at the output when the input displayed healthy unaltered activity (baseline, label = 0). A precision score of 1.0 means that every retrieved event is clinically-relevant. The algorithm for computing the precision score looped over the whole dataset. Hence the score took into account the spike emitted for all samples and by all four significant output channels. The ictal response latency was calculated as the time difference between the rising edge of the ictal event in the label signal and the rising edge of the first subsequent spike impulse in the recorded output. Specifically, the rising edge was the time stamp at which the signal crossed half of its maximum value. Same approach was adopted to compute the interictal response latency. Considering the first two spikes appearing at the output of the neuromorphic device after the onset of an ictal episode, their interspike interval ($ISI_{I-II}$) was estimated as the time elapsed between the falling edge of the first impulse and rising edge of the second one, procedure used also for each of the interictal discharges. Similarly, the interspike interval $ISI_{I-III}$ was computed as the time between the first and the third output spike, when applicable.

Figure 1 shows the results of a detection task performed by NET-TEN. Prerecorded LFP data of 100 s duration were converted into spikes in software using the SFE algorithm proposed by Kasabov et al [43] and delivered to the fabricated network through an arbitrary waveform generator. A second waveform generator synchronized with the first one was employed to carry the label signal and be able to display it on the oscilloscope together with NET-TEN neuronal spiking. Firing activity at the output signaled the detection of a pathological pattern in the LFP recordings.

Zoom In Zoom Out Reset image size

Figure 2 provides a summary of the designed hardware components. The neuromorphic system was implemented in a 180 nm CMOS technology node and its layout covered an area of 1.08 mm². The circuit voltage biases were adjusted by tuning the relative potentiometers on the PCB (depicted in figure 2(a)), while at the same time the corresponding firing activity generated at the output was compared with the label. In this way, it was possible to determine heuristically in which direction to steer each bias, in order to improve the classification accuracy. The fabricated chip is portrayed in figure 2(b). The voltage supply was set at 250 mV. The total static power consumption was 0.68 pW. The average power consumed by NET-TEN during a detection task calculated across samples was 3.50 nW. To evaluate the electrical behavior of the various functional modules that make up the neuromorphic network, these were fabricated as stand-alone components on a separate chip (figure 2(c)).

Zoom In Zoom Out Reset image size

Figure 3 illustrates the schematic of the three principal building blocks that constitute NET-TEN, namely the neuron, the excitatory post-synaptic current (EPSC) and the STDP circuits. As can be seen in figure 3(a), NET-TEN network has a feed-forward architecture composed by three sparsely-connected layers of ten neurons each. Every neuron of the input layer project to two neurons of the hidden layer, whose neurons also have a fan-in of two. Hidden and output layers are connected together in a one-to-one fashion, meaning every neuron of the hidden layer forms a synapse with one and only one neuron of the output.

Zoom In Zoom Out Reset image size

The neuron, depicted in figure 3(b), was initially proposed in [45] and further described in [37]. The circuit replicates biologically plausible dynamics, following the Izhikevich model. After performing the Monte Carlo analysis, the transistor sizes were modified with respect to the original implementation and the capacitance values were set to $C_\mathrm{v}$ = 76.86 fF and C_u = 614.88 fF, to alleviate the effect of process variations. One single neuron cell occupied an area of 1476.84 $\mu{m}^2$. In the schematic, the two state variables of the Izhikevich model, i.e. the transmembrane voltage $V_\mathrm{mem}$ and the slow recovery variable U, are represented by the voltage across the integrating capacitors $C_\mathrm{v}$ and C_u , respectively. $I_\mathrm{syn}$ is the synaptic current generated on the basis of the previous network layer and, when injected into the neuron, it perturbs its equilibrium. $M1-M3$ current mirror serves as a positive feedback to accelerate the rise of $V_\mathrm{mem}$ and bring it closer to the switching threshold of $M9-M10$ inverter. This increase is counteracted by the leakage current $I_{M4}$, proportional to U, which causes the probability of firing to be at a minimum immediately after a spike. Another current mirror, formed by the transistors $M1, M2$ and $M7$, couples U to $V_\mathrm{mem}$. A negative feedback introduced by transistor $M6$ discharges C_u , reducing U by an amount dependent on its own present value. The two inverters $M9-M10$ and $M11-M12$ function as a comparator with a fixed threshold, whose value derives from how the transistors are sized. When the transmembrane voltage exceeds the threshold, the comparator emits a negative pulse to turn on $M8$ and reset U, and a positive pulse to reset $V_\mathrm{mem}$ through $M5$. The reset mechanism is controlled by voltage biases $V_\mathrm{C}$ and $V_\mathrm{D}$, which are a means to alter the spiking pattern. Negative and positive pulses are harnessed at the network level to activate the synaptic circuits whenever a pre- or post-synaptic spike is elicited [37].

The EPSC circuit shown in figure 3(c) is responsible for inducing the post-synaptic current $I_\mathrm{syn}$ when triggered by a pre-synaptic spike and it constitutes the fundamental module of the silicon synapse. The EPSC is a low-pass filtered version of the original spike impulse, and presents a steep rising phase and an approximately exponential decay. In between spikes, the transistor $M_{E}3$ is on and conducts a current. Since at the same time $M_{E}1$ is off, no current can flow in this branch and therefore all $I_{M_{E}3}$ is diverted into $C_\mathrm{syn}$, precharging the node $V_\mathrm{syn}$ to $V_\mathrm{DD}$. When a pre-synaptic pulse arrives, the discharge phase starts. As $M_{E}1$ turns on and $M_{E}3$ deactivates, a current $I_\mathrm{w}$ is drawn from $C_\mathrm{syn}$. The amplitude of such current is directly proportional to the synaptic strength $V_\mathrm{w}$, so the greater is $V_\mathrm{w}$ the more $C_\mathrm{syn}$ is discharged. At the falling edge of the pulse pre, $M_{E}3$ begins to conduct current again, leading to a slow increase in $V_\mathrm{syn}$ value. The whole circuit operates in subthreshold, which means that the EPSC has an exponential dependence on $V_\mathrm{syn}$. The ratio between the sizing of $M_{E}1$ and that of $M_{E}3$ drives the time course of $I_\mathrm{syn}$, resulting in a fast onset followed by a slow decaying phase (corresponding to the recharge of $C_\mathrm{syn}$ after spike completion). The EPSC circuit response has been derived analytically in [37], which also illustrates the relationship between $I_\mathrm{syn}$ peak value and the synaptic weight $V_\mathrm{w}$. Simulations of the circuit with a power supply $V_\mathrm{DD} = 250$ mV, meaning the same one used in the experimental studies, and a static weight $V_\mathrm{w} = 200$ mV, which is the maximum bias value achievable with the tunable voltage dividers present on the PCB, led to an EPSC with peak current $I_\mathrm{syn_{peak}} = 486.39$ pA and a decay constant of τ = 212.83 $\mu s$. The layout of the EPSC was 462.09 $\mu{m}^2$ large.

The synaptic strength for each of the connections that link NET-TEN layers together is dynamically updated in agreement with the STDP rule and it can take up any analog value between $V_\mathrm{DD}$ and ground. NET-TEN featured the STDP CMOS implementation conceived by Indiveri [44], which was here redesigned to work at ultra-low power (figure 3(d)). Amplitude and temporal evolution of potentiation and depression could be regulated respectively by voltage biases $V_\mathrm{tp}$ and $V_\mathrm{td}$, and $V_\mathrm{pb}$ and $V_\mathrm{db}$. The new implementation occupied 1054.75 $\mu \mathrm{m}^2$. Given that the synaptic weight was stored on a capacitor $C_\mathrm{w}$ = 206.30 fF, the synapses were able to retain the computed value only for a short period of time, in the order of tens of milliseconds.

The analyzed dataset comprises MEA recordings from the parahippocampal cortex (CTX) of n = 10 brain slices from n = 10 mice, to ensure fair representation of biological variability. The extracted samples were 100 s long signal segments around an ictal event. Three samples also featured interictal discharges, while one of them only contained an artifact. The dataset content is summarized in table 1. Figure 4 depicts a hippocampus-CTX slice coupled to a MEA and representative of the epileptiform patterns recorded by the electrode marked by the blue circle. During the measurements, the UP spikes were fed to five input neurons, the DW spikes to the other five, as shown in figure 3(d). All ten output pins were probed with the oscilloscope, but only four out of these ten were deemed to be informative for pattern recognition. In fact, despite all neuron and synapse circuits share the same nominal biases values, device mismatch affects the individual instances behavior. The described results refer to said four output channels.

Zoom In Zoom Out Reset image size

The output response exhibited by NET-TEN in three distinct cases can be observed in figure 5. The precision score was 0.998. The device fired in correspondence of ictal events with a mean latency of 110.4 ms, while it took about 126.6 ms to respond to an interictal discharge. Table 1 details the response latency calculated for each LFP sample and for each retrieved output signal. Observe how the latency varied among the output channels, with some channels more reactive and others less reactive to the stimulus, depending on the internal dynamics of the network. Further classification capabilities were explored, looking for elements that would allow a distinction between ictal and interictal occurrence, and possible ways to discard artifacts. The amount of spikes released at the output differed between classes, and so did the interspike interval ($ISI_{I-II}$) in-between the first two spikes. The median spike count for ictal events amounted to 1117.5 with an interquartile range (IQR) of 1762.5. For interictal activity, the median spike count was 4.5 with IQR 7. In the case of the artifact, median and IQR were 27.5 and 22.8, with the first quartile being 19.3 and the third quartile equal to 42. Coming to the $ISI_{I-II}$, the values were first averaged across output channels. The median $ISI_{I-II}$ was 4.3 (3.7–4.8) ms for ictal and 7.8 (3.8–33.4) ms for interictal patterns. Instead, since there was only one sample including one artifact, the outcome of the average of $ISI_{I-II}$ across measured output signals came down to a single value corresponding to $ISI_\mathrm{artifact} = 562\,ms$. Figure 6(a) provides a visual representation of the feature space, where the number of output spikes is on the X-axis and $ISI_{I-II}$ on the Y-axis. In this plot both measurements are averaged across output channels, so each data point corresponds to an occurrence, either an ictal event, an interictal discharge or the artifact. As you can see, the three classes are clustered in three distinct groups. In particular, the clinically relevant events have a lower $ISI_{I-II}$ compared to the artifact, which means their onset is associated with a higher instantaneous firing frequency at the output. However, plotting the measurements relative to the different channels separately, as in figure 6(b), suggests that the use of $ISI_{I-II}$ could lead to a misclassification of the artifact. The interval $ISI_{I-III}$ is a better means of discrimination (figure 6(c)). Also, to take full advantage of the redundancy given by the multiple output channels, a potential postprocessing algorithm could look for the minimum $ISI_{I-III}$ across channels (see figure 6(d)). Interictal and ictal instances can be further discerned based on their spike count.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Listing 1 shows the pseudocode of the algorithm used to isolate the ictal events in NET-TEN's output signal. For each pulse in the signal, the algorithm calculates the interval between its falling edge and the rising edge two pulses ahead, and then compares it to a reference value set at $ISI_{thr} = 10$ ms. The reference value was extracted from figure 6(d). If the condition $ISI_{I-III} \lt ISI_{thr}$ is verified, which rules out the artifact, the algorithm checks whether 10 more spikes occur in the following 190 ms. If this second condition is also satisfied, the instance is tagged as an ictal event. The detection mark is placed in correspondence of the last spike of the sequence (the 13th spike). The ictal event is considered detected when it is tagged in at least one output channel. Table 2 lists the detection delays obtained for the different samples, quantified as the time between the detection mark and the ictal onset in the relative label signal.

Zoom In Zoom Out Reset image size

An almost equivalent way to implement the above-mentioned two-step algorithm online and on-chip is by using a counter to count the number of spikes in the two predetermined time windows. A 8-bit counter, implemented in a 180 nm CMOS technology node, with clock frequency of 1 kHz can serve to continuously count up to 10 ms and 190 ms. According to the simulations, the dynamic power consumption of such counter amounts to 31 pW at 300 mV voltage supply. Meanwhile, a 3-bit counter counts the spikes fired at each output channel. The power consumption of this counter is signal-dependent, so its estimation is not straightforward. In terms of energy, the counter consumes 11.4 fJ/spike, which translates into 11.4 fJ/spike × 13 spike/classification = 148.2 fJ/classification.

Table 3 showcases the specifications and the performance of the developed device compared with other state-of-the-art systems for seizure detection.

Table 3. Comparison with state-of-the-art miniaturized systems for seizure detection.

	Salam et al TBCAS 2011 [16]	Yoo et al JSSC 2012 [46]	Chen et al JSSC 2014 [19]	Altaf et al TBCAS 2016 [18]	Sharifshazileh et al Nat. Commun 2021 [35]	This work
Recorded signal	iEEG	scalp EEG	iEEG	scalp EEG	iEEG (HFO)	LFP
Technology	180 nm	180 nm	180 nm	180 nm	180 nm	180 nm
Area	2 mm2 $\dagger$	$5 \times 5$ mm2 $\dagger$ (8 channels)	$2.76 \times 4.88$ mm2 $\dagger$ (8 channels + stimulation)	$5 \times 5$ mm2 $\dagger$ (8 channels)	77.2 mm2	1.08 mm2 $\ddagger$
Supply	1.8 V	1.0 V (LSVM)	1.8 V	1.0 V (NLSVM)	1.8 V	0.25 V (NET-TEN) 0.3 V (counters)
On-chip classification	Voltage level and high frequency detectors	Linear support vector machine	Bio-signal processor (Entropy and Spectrum)	Non-linear support vector machine	Spiking neural network	Spiking neural network
Power consumption	44.85 nW (detector)	66 µW + 2.03 µJ/classification	162.31 µW/channel (BSP)	100.6 µW (NLSVM) + 1.83 µJ/classification	555.6 µW	3.53 nW + 148.2 fJ/classification
Detection delay	13.5 s	${\lt}2$ s	0.8 s	2 s	${\sim}30$ ms	56.7 ms [79.1 ms] (NET-TEN) 0.23 s [2.13 s] (ictal only)

Includes analog front end.Does not account for counters of postprocessing algorithm.HFO = high frequency oscillations, BSP = bio-signal processor.LSVM = linear support vector machine, NLSVM = non-linear support vector machine.