Phase Synchronization Operator for On-Chip Brain Functional Connectivity Computation

IEEE TRANSACTIONS ON BIOMEDICAL CIRCUITS AND SYSTEMS, VOL. 13, NO. 5, OCTOBER 2019 957 Phase Synchronization Operator for On-Chip Brain Functional Connectivity Computation Manuel Delgado-Restituto , Senior Member, IEEE, James Brian Romaine , and Ángel Rodríguez-Vázquez , Fellow, IEEE Abstract—Thispaperpresentsaninteger-baseddigitalprocessor for the calculation of phase synchronization between two neural signals. It is based on the measurement of time periods between two consecutive minima. The simplicity of the approach allows for the use of elementary digital blocks, such as registers, counters, and adders. The processor, fabricated in a 0.18-µm CMOS process, only occupies 0.05 mm2and consumes 15 nW from a 0.5 V supply voltage at a signal input rate of 1024 S/s. These low-area and low-power features make the proposed processor a valuable computing element in closed-loop neural prosthesis for the treatment of neural disorders, such as epilepsy, or for assessing the patterns of correlated activity in neural assemblies through the evaluation of functional connectivity maps. Index Terms—Functional connectivity, low power CMOS VLSI, neural signal processing, phase synchronization, seizure detection. I. INTRODUCTION OVER the past decade, there has been a growing interest in the quantification of functional connectivity between measurements of neural activity. This refers to the statistical evaluation of coupling strengths between brain regions to identifythose involvedin specifictasks,includingsensoryresponses, motor responses and intellectual or emotional processing. Functional connectivity can be assessed on neural signals captured at different spatial resolutions, from single-unit responses acquired by micro-electrode arrays to aggregated signals obtained through electroencephalography (EEG) or functional magnetic resonance imaging (fMRI) measurements, and it can be calculated on the basis of linear (e.g., correlations and coherence measures [1]), nonlinear (e.g., generalized synchronization [2]) or information-related techniques (e.g., cross mutual information [3])—see [4] for an in-depth review. Functional connectivity does not determine the specific direction of information flow but only shows the regions which are most likely connected. To explicitly assess the influence that one region exerts over Manuscript received April 6, 2019; revised May 14, 2019 and June 24, 2019; accepted July 19, 2019. Date of publication July 29, 2019; date of current version November 4, 2019. This work was supported in part by the Spanish Ministry of Economy and Competitiveness under Grant TEC2016-80923-P and in part by the ONR under Grant ONR N00014-19-1-2156 and the FEDER Program. This paper was recommended by Associate Editor A. Basu. (All authors contributed equally to this work.) (Corresponding author: James Brian Romaine.) The authors are with the Institute of Microelectronics of Sevilla, University of Seville-CSIC, Sevilla 41092, Spain (e-mail: [email protected]; [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TBCAS.2019.2931799 another, i.e., to derive the effective connectivity [5] and, hence, to reveal the subjacent drive-response map of the neural system, more complex techniques based on the analysis of long time series are typically used. They include the Granger-causality [6], Directed Transfer Function (DTF) [7] or Partial Directed Coherence (PDC) [8]. Besides providing means to expose communication networks in the brain, clinical studies have demonstrated that functional connectivity also gives relevant information to distinguish between normal and pathological brain states, in particular, when connectivity is estimated by synchrony evaluation techniques [9]–[11]. Indeed, it has been shown that abnormal synchronizationpatternsareassociatedwithdifferentneurologicaldisorders, such as, epilepsy [12], Alzheimer’s disease [13], Parkinsons disease [10] or schizophrenia [14]. Moreover, changes in functional connectivity, possibly together with other univariate measures extracted from neural recordings, have been found useful for the early detection of epileptic seizures with good sensitivity and specificity [2], [15], [16]. This later feature is actually exploited in recent closed-loop neural prostheses for treating epilepsy [17]. These prostheses detect specific bio-markers in neural signals during the pre-ictal stages of seizures and, following detection, regions of the brain are stimulated for preventing the spreading of seizure activity. For instance, the observed desynchronization at the initiation of seizures in human intracranial and extracranial recordings [18] has been used as a bio-marker for predicting ictal periods long before a seizure actually occurs [19]. The referred applications on neurophysiology, disease monitoring and therapy suggest the development of dedicated integrated circuits for functional connectivity evaluation, in particular, in those scenarios demanding power efficient solutions (e.g., implantable neural prostheses) or fast on-site operation to avoid neural activity transfer and off-line computations by a host computer (e.g., neurostimulators for seizure abortion). However, not all the candidate algorithms for derivingfunctional connectivity lead to hardware friendly integrations. Most often, they require the storage of long sequences often combined with complex operators, e.g., Fourier Transform modules in the case of coherence measures [1] or entropy calculators in the case of information-related techniques [3]. Hence, their integration, even assuming extensive usage of parallelism for addressing the connectivity between multiple channels, would demand large area occupation. Similarly, synchrony evaluation often relies on analytical techniques involving trigonometric functions, which 1932-4545 © 2019 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. 958 IEEE TRANSACTIONS ON BIOMEDICAL CIRCUITS AND SYSTEMS, VOL. 13, NO. 5, OCTOBER 2019 demandbulkyandpower-hungrycomputations[20]whichmake multichannel extension difficult as scalability is reduced. This paper addresses the integrability of functional connectivity operators in silicon and presents a simple, hardware friendly algorithm and a companion CMOS implementation as a demonstrator. The algorithm focuses on the measurement of synchrony between two signals and, similar to [21], relies on the identificationof events withno need to evaluatethe instantaneous phaseof thesignals.Althoughthealgorithmcanbeextendedtoanyneural recording and identifiable event (e.g., action potentials in single unit measurements), herein it is illustrated with narrow-band filtered surface or intracranial EEG signals, using local minima asevents[22],[23]. The demonstrator uses integerarithmeticfor calculations and consists of counters, registers and comparators, as main building blocks. In spite of its simplicity, the proposed operator has an accuracy comparable to that of more complex algorithms, such as the the mean phase coherence [24], but with much lower area and power consumptions. Indeed, the processor, fabricated in a 0.18 μm CMOS technology, consumes 15 nW from 0.5 V supply and only occupies 0.05 mm2. The paper is organized as follows. Section II defines the phase synchronization concept and Section III briefly describes some of the methods most commonly used for the on-chip quantification of synchrony between two signals. Then, Section IV introduces the proposed algorithm and Section V presents the structure and operation of the developed processor. Section VI describes the methods followed for testing and shows the achieved performance. It presents experimental results, including functional connectivity maps, obtained by driving the integrated operator with neural recording streams available in an epilepsy database [25]. Comparisons with simulations of analytical techniques for synchrony computation and seizure detectionarealsogiventoillustratetheaccuracyoftheapproach. Finally, Section VII presents some concluding remarks and future plans. II. PHASE SYNCHRONIZATION CONCEPT Phase Synchronization (PS) is expressed by the phase locking condition n·ϕs1−m·ϕs2=const., where nand mare integers, and ϕs1and ϕs2denote the instantaneous phases of the oscillating signals s1(t)and s2(t). In the presence of noise, this condition is relaxed to the weaker constraint |n·ϕs1−m·ϕs2|<const.or to the frequency entrainment condition, n·ωs1−m·ωs2=0, where ωs1and ωs2denote the angular frequencies of s1(t)and s2(t), respectively, and · denotes averaging over time [26]. The measurement of PS over time typically encompasses two main steps: first, the estimation of phase angles ϕs1and ϕs2and, second, the quantification of the amount of synchronization over a given period of time. This period is usually defined by a sliding window of Nsamples and, hence, the quantification inherently carries out a data reduction process. Synchronization measures often contain rapid fluctuations. Hence, an additional smoothing processcanprovideamoreusefulbaseline for assessingchanges in synchrony. In practical situations, PS estimation is preceded by the pre-filtering of the neural signals such that s1(t)and s2(t)are band-limited into specific frequency bands. For instance, in the case of signals obtained by EEG or electrocorticography (ECoG)methods, these bandsare conventionallydefinedas: δ∼ 0.1−3Hz,θ∼4−7Hz,α∼8−12 Hz, β∼12−30 Hz, γ∼ 30−60 Hz and high-γ>60 Hz. Given the narrow-band constraint and considering that signals are obtained from the same physiological system (i.e. the brain) with similar recording techniques (i.e. surface or intracranial EEG), it will be assumed in the following that PS can be more likely found for phase locking integers n=m=1[24], [26]. III. REVIEW OF ON-CHIP PS SOLUTIONS One popular approach for the calculation of the instantaneous phase of a signal relies on the construction of the analytic representation sa(t)=s(t)+j˜s(t), where ˜s(t)is the Hilbert Transform of s(t)[see (A.1) in the Appendix]. Using time-frequency domain techniques, the Hilbert Transform can be derived by [31] (i) obtaining the Fourier transform of s(t), (ii) nulling the negative frequency components, (iii) calculating the inverse Fourier transform and(iv) taking the imaginary part of the result. Given the non-stationary nature of neural signals, the Fourier transform should be applied over short intervals. This suggests the use of the Short Time Fourier Transform (STFT) algorithm instead of the classical Fast Fourier Transform (FFT) [4]. A simpler method for the generation of the analytic signal representation sa(t)uses Hilbert Transformers based on digital filtering. These transformers can be implemented in complex domain or formed by the combination of two sub-components; a Hilbert filter with an ideal transfer function −jsgn(f)for the generation of the imaginary part of sa(t), and a delay equal to that of the Hilbert filter for the real part. In both cases, the total group delay of the signal components should be kept constant as to not cause phase errors. The vector analyzer presented in [27] follows the second implementation approach and employs 16-tap Hilbert finite impulse response (FIR) filters for obtaining ˜s(t), and 16-tap All-Pass FIR filters for implementing the delay. In a more recent implementation [28], the in-phase and quadrature components are directly provided by the mixed-signal acquisition front-end, and four 64-tap FIR filters are used for band selection. Further, [34] introduces an improved implementation with 9×area reduction by means of resource re-utilization. Another transform-based approach for deriving the instantaneous phase of a signal uses Wavelets [see A.4 in the Appendix]. In this case, the on-chip implementation of the method requires Multiply-and-ACcumulate (MAC) processors to correlate the input time series with wavelet templates. As an example, [29] proposes a charge-recycling mixed-signal MAC processor for epilepsy detection in which time series of 1024 samples are serially loaded and correlated in parallel with all the Morlet wavelet templates (coded in 32 frequency bins with 4-bit coefficient resolution) stored in the on-chip memory. The trigonometric functions involved in (A.2) and (A.3) as well as the computation of the phase locking value in (A.5) can be readily implemented by CORDIC (COordinate Rotation DIgital Computer) units. This is actually done in [20] where a digital signal processing (DSP) block formed by three 16-bit CORDIC cores and two 32-tap moving average FIR filters are DELGADO-RESTITUTO et al.: PHASE SYNCHRONIZATION OPERATOR FOR ON-CHIP BRAIN FUNCTIONAL CONNECTIVITY COMPUTATION 959 TABLE I PERFORMANCE SUMMARY OF PROPOSALS RELATED TO ON-CHIP PS CALCULATION aPower consumption mainly due to the 4 ×128 8-bit ΣΔ algorithmic ADCs. At 15 kHz parallel sample rate, the bank of ADCs dissipates 6.3 mW from a 3.3 V supply [30]. used for computing the phase locking value (PLV) from the analytical signal sa(t). In [35], FIR filters are replaced with an IIR approximation resulting in a 60% decrease in group delay latency. Another alternative approach for computing PLV in which trigonometric functions are replaced by linear algebra calculations is presented in [36]. Numerical methods for computing PS, such as the one based on probabilities of recurrence, pr(τ), described in the Appendix, typically use simpler logic than transform-based techniques at thepriceof increased latencyand morememoryresources.Thus, for instance, the evaluation of pr(τ)in (A.7) requires the collection and combination of phase space trajectories, demanding for a large memory allocation. Nevertheless, it is worth observing that,depending onτ, thereare coincident terms inthe sum(A.7). This opens doors for reducing the computational complexity of the algorithm by avoiding recalculations. This is actually the approach followed in [32], where a differential scheduling for the computation of the correlation integral of a signal (similar to the calculation of recurrences) allows for a substantial reduction in complexity. Following the calculation of pr(τ), a MAC unit should be used for the calculation of the cross-correlation in (A.8) to serve as a synchronization index. As shown in (A.9) and (A.10), numerical methods can also be used for the computation of synchronization indexes. Although the computational complexity is relaxed as compared to the PLV (in particular for the index based on Shannon entropy), long sliding windows are needed to populate the histograms and obtain a meaningful distribution of probabilities. Hence, circuit simplicity is counterbalanced with an increase in latency and memory depth. Table I illustrates the performance of the different integrated circuits presented in this section, together with the performance of the presented work (to be discussed in Section VI). The two first rows represent for complete PS processors while the rest corresponds to building elements which can be eventually used depending on the particular PS implementation strategy. The table clearly shows that the algorithms used for PS calculation, although easily programmable in computers, tend to occupy large area and consume significant power when dealing with silicon integration. The problem is even more acute in the context of functional brain connectivity, where multiple phase synchronization indexes from different sites are needed. For instance, if functional connectivity has to be measured between 16 recording sites, as it is typically done in a standard 10–20 EEG electrode system, 120 different combinations of neural signal pairs should be simultaneously addressed. This means application specific integrated circuits (ASICs) for functional connectivity quickly become far to big and/or power hungry, even if sophisticated multi-thread sharing DSP techniques are employed. IV. EVENT-BASED PHASE SYNCHRONIZATION In this section, an algorithm for phase synchronization computation which aims to alleviate the area/power consumption burden for silicon integration of previous approaches is proposed. Similar to the transform-based methods, the proposal is also based on the estimation of instantaneous phases in band-limited signals, however, instead of using complex filtering or wavelet convolutions, phases are extracted through the identification of distinct marker events in waveforms. Without loss of generality, in this work it is assumed that such events correspond to the local minima of the signal at time instants tn. The time interval between two consecutive minima correspond to one complete cycle and, therefore, the phase increment during this period is exactly 2π. Hence, the signal phase can be approximated for any arbitrary time instant tn<t<t n+1 as, ϕs(t)≈2πt−tn tn+1 −tn +2πn (1) where the term 2πn unwraps the phase angle and Ts n=tn+1 − tnmeasures the duration of the transition period between the minima, as illustrated in Fig. 1. Note that this approximation notably simplifies phase estimation as the problem basically reduces to measure time intervals. For this reason, the proposed algorithm is denoted as Delay Difference Analysis (DDA). Based on the approximation in (1), an index for quantifying thesynchronizationlevelbetweentwosignalshasbeenproposed in [21]. It relies on creating histograms on how many times the events in one signal are preceded by events in the other. The approach has low computational cost however, similarly to the numerical techniques described in the Section III, it needs longtimeseriestobuildmeaningful statistics and, in somesense, it wastes the time information between events. The proposed approach for measuring synchronization is similar to the PLV technique in (A.5) but avoids the use of trigonometric functions. DDA analyzes the time difference between the transition periods of both signals, not their relative 960 IEEE TRANSACTIONS ON BIOMEDICAL CIRCUITS AND SYSTEMS, VOL. 13, NO. 5, OCTOBER 2019 Fig. 1. Transition periods and approximated phases in two signals. Starting at point A, the DDA algorithm generates the following absolute differences between transition periods: |T1 1−T2 1|,|T1 2−T2 2|,|T1 3−T2 3|and |T1 5−T2 4|. phase ϕ, and, remarkably, it operates on-the-fly with no need to store long signal sequences. The procedure simply consists in pairing transition periods from both signals as soon as they are detected. If a pair is formed between the i-th transition period of signal s1(t)and the j-th transition period of s2(t), the absolute difference Tn=|T1 i−T2 j|is calculated, where T1 iand T2 j denote, respectively, the durations of said transition periods. If two or more transition periods in one signal are detected before a transition period is found in the other, only the last transition period of the former signal is accounted for in the computation of Tn. Of course, this situation is more or less likely to happen depending on the band limitations of the input signals. For example, in the theta band which has a narrow frequency spread, the probability of multiple transitional detection is less than that of the beta band. Fig. 1 shows an example of the procedure. Let us denote by Kthe number of absolute differences Tn computed in the k-th sliding window of length N, and let us assume with no loss of generality that N=2 m. A synchronization index between signals s1(t)and s2(t)can be defined as, S(r) dda(k)=1−2r Nmin K−1  n=0 |Tn|−Tos,N 2r(2) where r∈[0,m]is a user-defined selectivity control parameter and Tos is a positive integer offset for adjusting the synchronization index range. Fig. 2(a) plots S(r) dda for Tos =0and different selectivity configurations assuming that s1(t)is a sine wave at 20 Hz,ands2(t)isachirpwhosefrequencylinearlysweepsfrom 10 to 30 Hz. Both signals are sampled at 1024 S/s and quantized at 10-b resolution. Every calculated index is obtained assuming a non-overlapping sliding window of length N= 1024. Clearly, the higher the selectivity parameter, the narrower the band for which the synchronization index obtains non-null values. No matter the rparameter, the plots in Fig. 2 are not symmetrical around the tone frequency of s1(t); they decay linearly for chirp frequencies below said tone but fall away more smoothly Fig. 2. Synchronization indexes obtained through the DDA algorithm for different selectivity parameters, using a sine wave and a swept-frequency cosine as input signals. (a) Tos =0;(b)Tos =6. for higher frequencies. This is because the number of absolute differences Kis limited by the smaller frequency, fmin, of both signals as K=N·fmin/fs=Tout ·fmin, where Tout is the throughput period of the DDA index. In Fig. 2(a), it is also observed that the peak value of the synchronization index decreases with the selectivity parameter. This can be compensated, without drastically increasing the sampling frequency or the number of samples per observation window, by means of the offset parameter Tos. This is illustrated in Fig. 2(b) where Tos =6 making the synchronized indexes for all selectivity parameters mostly cover the range between 0, indicating no coupling, and an approximate 1 for perfect coupling. Note that too large of an offset may eventually make the synchronization index larger than 1, hence, a limiter, which can be easily implemented by an overflow detector, should be used together with (2). For comparison purposes, Fig. 3 shows the PLV index (dashed line) for the same experiment as in Fig. 2 together with the index S(4) dda with no offset. Both approaches give comparable results,withtheDDAproducingapiecewiselinearapproximate of the main lobe of the PLV index and eliminating the smaller side lobes. This is reasonable because the PLV index can be approximated as PLV ≃1+αN−1 k=0 |ϕk|/N for small ϕ values,whereαisafittingparameter(see(A.5)intheAppendix). Given the relationship between phases and time delays as expressed in (1), a clear connection between the phase locking value and the Delay Difference Analysis can be established. DELGADO-RESTITUTO et al.: PHASE SYNCHRONIZATION OPERATOR FOR ON-CHIP BRAIN FUNCTIONAL CONNECTIVITY COMPUTATION 961 Fig. 3. Results from DDA and Hilbert transform approaches for the same experiment in Fig. 2. The inset shows a zoom of the shaded area. Fig. 4. Block diagram of the DDA circuit. V. DDA CIRCUIT DESIGN Fig. 4 shows the block diagram of the circuit implementing the DDA algorithm. Initially both signals (assumed bandpass limited off-chip) are smoothed by simple exponential filters (to be discussed later), and, afterwards, they are passed through corresponding event detectors for identifying the time instances leading to the calculation of the transition periods. As mentioned, local minima will be considered as marker events in this work. As shown in Fig. 4, each detector consists of (i) a digital Fig. 5. (a) Combinational circuit for the minima logic in the event detection blocks of Fig. 4. (b), Example of neural waveform and the detected minima (encircled) based on the event detection system. comparator driven by the current and previous samples of the signal, sx(i)and sx(i−1), (ii) a shift register, which stores the Mmost recent results of the pair-wise comparison and (iii) a combinational block which, based on the sequence stored in the register, decides if a minimum is present. This latter block, shown in Fig. 5(a), allows for a finite number Qof potential outliers in the up and down transitions of the signal. This is illustrated in Fig. 5(b). In spite of a short transition upwards in the second event, the detection system is still able to identify the correct local minimum. In the proposed implementation M=10and Q=2. Based on the detected events, a finite state machine calculates the differences between transition periods, Tn.Asshownin Fig. 4, the finite state machine uses up-counters to determine the number of clock cycles between two consecutive minima in a signal. When a pair of minima are detected, the output of the corresponding counter is stored in latches and, afterwards, the counter is reset to start a new count. If a transition period in s1(t) is preceded by a transition period in s2(t), the synchronization indexing block calculates the absolute difference between the two.Otherwise,ifthetransitionperiodins2(t)isyettobestored, thefinitestatemachinecontinuessearchingfortransitionperiods in s1(t), overwriting the values in the latches until a transition period in s2(t)is detected. The same procedure holds if the role of s1(t)and s2(t)are swapped. Anaccumulatorinthesynchronizationindexingblockobtains the sum of the absolute differences between transition periods until the number of samples per observation window, N,is reached. At this point, the synchronization index is obtained according to (2). In practice, to avoid floating point arithmetic, a scaled version of the index N·S(r) dda(k)is used. Hence, only comparators, shift registers and digital adders are needed for the calculation. Similar to that of the input signals at the pre-processing stage, the synchronization index is smoothed by a simple exponential 962 IEEE TRANSACTIONS ON BIOMEDICAL CIRCUITS AND SYSTEMS, VOL. 13, NO. 5, OCTOBER 2019 Fig. 6. Block diagram of the smoothing block. Fig. 7. Correlation of the DDA synchronization index and the PLV processor in [20] with an ideal noise-free PLV calculation when the input signals are contaminated by Gaussian noise. filter for removing fast fluctuations [23]. This filter is defined by the recursive equation: y(k)=(1−α)·y(k−1) + α·S(r) dda(k)(3) where y(k)is the filtered output and αis a smoothing factor comprised between 0 (maximum smoothing) and 1 (no smoothing). In practice, this parameter is chosen as a negative power of two, so that multiplications reduce to simple shifting operations in a register. Fig. 6, shows a block diagram of the filter which is composed of only two shift right registers, which are responsible for the power of two smoothing factor, a full subtracter and a full adder.A similar smoothingapproachhas been recentlyproposed in [35]. Fig. 7 illustrates the noise performance of the DDA circuit. It shows the correlation between the ideal PLV values obtained with two noise-free 60 minutes long neural signals, and the synchronization indexes, derived from the DDA approach, when the same signals are contaminated with white Gaussian noise. The correlation is plotted against the spot Signal-to-Noise Ratio (SNR), assuming 1 Hz bandwidth, at the maximum frequency of the βband, i.e. 30 Hz. PS estimations are also derived on the same frequency band. A Matlab model, closely matching the circuit architecture in Fig. 4 for N= 1024,M=10,Q=2 and α=1/32, has been used in the analysis, assuming 10-b resolution and a sampling rate of 1024 S/s. Parameters rand Tos have been adjusted for the best fit to the ideal PLV response. Spot SNRs in the range between 10 dB and 70 dB at 5 dB steps have been considered and, for each SNR, a Montecarlo analysis (100 instances) has been carried out. The error bars show both the mean value and the ±3σof the synchronization indexes. It is worth observing that for spot SNRs above 30 dB the average correlation between the DDA approach and the noise-free PLV technique is larger than 90%. Fig. 7 also shows the correlation between the noise-free PLV values and those obtained with a 16-tap Hilbert transformer with passband ripple of 0.01 dB when input signals are corrupted by noise, following the implementation in [20]. Similar to that of the DDA case, performance worsens for SNRs below 20 dB, although the DDA obtains better results with very noisy recordings thanks to the outlier suppression technique used in the minima detectors. Table II shows the mean value and the standard deviation of the correlation between the PS values obtained with the DDA processor with different parameter configurations and the ideal PLV routine for a spot noise of 30 dB. Note that the correlation performance is particularly sensitive to the smoothing factor α and that only modest improvement is observed by increasing the resolution of the input signals. Parameters Mand Qmainly affect PS calculations at low SNRs. The configuration marked in bold has been selected for implementation as it offers a good trade-off between hardware complexity and performance. To assess the impact of the noise behavior shown in Fig. 7, it must be emphasized that surface and intracranial EEG signals exhibit 1/fx-like power spectra at low frequencies [37], [38]. This makes low frequency bands more tolerant to the aggregated noise contributed by biological tissue, electrodes and recording front-ends [39], [40]. Indeed, it has been shown that even with modest amplifier noise specifications, without suppressing flicker contributions, SNRs larger than 30 dB can be obtained at low frequencies [41]. This suggests that PS calculations can be relevant as a neurological biomarker in the βand lower frequency bands although the reliability decreases for the high-γband and higher frequencies. This observation is indeedconfirmedinSectionVI,whereit is shownthat DDAand PLV obtain similar PS results, while the DDA approach offers light and modular hardware implementations, easily scalable to large neural recording arrays. VI. EXPERIMENTAL MEASUREMENTS Fig. 8 shows a microphotograph of the chip implementing the DDA algorithm. The chip has been fabricated in a 0.18 μm CMOS process and occupies an active area of 0.05 mm2in- cluding the main processor core, two series-to-parallel input (SPI) ports, one parallel-to-series output (PSO) port and some additional buffers for testing intermediate nodes. The I/O ports work 10 times faster than the rest of the ASIC in agreement with the 10-b resolution used for the input signal samples. The total number of integrated gates is 6053. Before integration, the processor core was verified on an ARTIX 7 FPGA. Only 92 slices were required for implementing the proposed algorithm. DELGADO-RESTITUTO et al.: PHASE SYNCHRONIZATION OPERATOR FOR ON-CHIP BRAIN FUNCTIONAL CONNECTIVITY COMPUTATION 963 TABLE II DDA CORRELATION TO IDEAL PLV MEASURED AT 30 DBSPOT SNR FOR DIFFERENT CIRCUIT CONFIGURATIONS Fig. 8. Circuit layout fabricated in a 0.18 μm CMOS process and occupies an active area of 0.05 mm2. Fig. 9. Laboratory functional test setup, showing the logic analyzer, power supply and the laptop used to transfer data files. The ASIC uses circuit elements from a full-custom digital library for 0.5 V power supply operation, based on low-leakage transistors in weak inversion. The library blocks were generally designed using conventional CMOS logic, with device dimensions close to minimum size, as they were found to offer an acceptable trade-off between power dissipation and operating frequency [42]. Indeed, thelibrary was exhaustively verified and characterized under PVT deviations assuming a maximum clock frequency of 25 kHz, i.e., about twice the operation frequency needed by the I/O ports of the ASIC. In spite of the reduced noise margin and the slow switching transitions, post-layout simulations of the chip showed complete functionality with no data loss for 0.5 V operation. Fig. 9 shows the test setup used for functional verification. The chip was mounted on a dedicated PCB together with a set Fig. 10. Illustration of the system response versus supply voltage for two input tones at 25 and 30 Hz. of adjustable level shifters and regulators for scaling the logic levels and supply voltage of the ASIC. Low voltage operation was experimentally verified by driving the chip with different tones in the βband and comparing the generated synchronization index stream with electrical and behavioral simulations. A logic analyzer (Agilent 16823A) was used for synthesizing the tones, sampled at 1024 S/s, and for serially retrieving binary information from the output and other test points of the ASIC through a pod connector. For two tones at 25 and 30 Hz (100 k samples each), deviations from the ideal synchronization index were noticeable for supply voltages below 0.35 V. At this point, the index starts to drop and, for 0.2 V supply, the system breaks down completely with no observed activity. The reset signal generatedintheASICafterand-ingthelatchedoutputsofthetwo minima logic blocks (see Fig. 4) was also monitored along this transition. This signal goes high every time a transition period is computed. For supply voltages above 0.35 V, the train generated by the ASIC is quite regular (only minor variations due to the quantization of the inputs) in good agreement with the expected response, however, pulses are eventually missed below 0.35 V and, at 0.2 V supply, they are no longer observable. Electrical simulations show that this behavior is compatible with data transfer failures from the SPI input ports of the ASIC. This is illustrated in Fig. 10 which plots: (i) the experimentally observed degradation of the synchronization index as the supply voltagedecreasesaccordingto Smeas/Sid,whereSmeas andSid denote, respectively, the measured and ideal values of the index and, (ii) the simulated SPI register position, normalized to the vector resolution, at which input data fails to be shifted. In both cases, 1 indicates no error while 0 means complete breakdown. It is worth observing the good agreement between both curves. Similar behavior was obtained for other tones in the βband. At a signal input rate of 1024 S/s, a throughput period of 1 s and 0.5 V supply voltage, the power consumption of the chip core is 15 nW. Table III shows the power budget of the chip. Event detection and time stamp calculations by the finite state machine are the two most power demanding tasks of the processor due to the higher switching activity. Table IV shows the normalized performance of the ASIC compared to other designs in Table I. Variables have been normalized to a 964 IEEE TRANSACTIONS ON BIOMEDICAL CIRCUITS AND SYSTEMS, VOL. 13, NO. 5, OCTOBER 2019 TABLE III POWER BREAKDOWN OF THE SYNCHRONIZATION PROCESSOR BASED ON THE BLOCKS IN FIG.4 TABLE IV NORMALIZED PERFORMANCE SUMMARY Fig. 11. Test setup showing the ASIC mounted on a dedicated PCB, mbed LPC1768 micro controller and dedicated level shifters. single channel and, in the case of core area occupation, an additional normalization has been applied by arbitrarily assuming all chips are fabricated in CMOS 65 nm. Accordingly, the original area per channel has been scaled by (65 nm/Wmin)2, where Wmin is the feature size of the fabrication technology. Note that the proposed implementation exhibits the smallest normalized area and it is among the designs with lower energy consumption. A. Experimental Setup With Neural Recordings The processor has also been verified and validated using neural recording data available in the European Epilepsy database [25]. This database contains recordings of over 275 patients including 225 scalp recordings and 50 intracranial recordings. As very long recording sessions, often lasting for more than 6 h, are included in the database, a specific experimental set-up has been devised. As shown in Fig. 11, it uses an mbed LPC1768 MCU based on a 32-bit ARM Cortex-M3 core from NXP Semiconductors. The test platform has been expanded with an application board which includes an RJ45 Ethernet connector and a USB-A Host/Device port (in addition to the built-in USB drag’n’drop FLASH programmer). The MCU defines settings and transmits driving stimuli to the DDA chip through the available GPIO ports and other I/O interfaces. Neural input signals are stored as CSV files in a 32 GB Fat-32 formatted stick memory plugged into the USB-A port of the test platform. These files are sequentially read by the MCU and serially transmitted to the SPI ports of the ASIC. Long testing sessions have been run uninterruptedlyusingthissetupwithnoneedtosplitthedatabase recordings into smaller blocks (this would clear the internal states of the ASIC). The synchronization index generated by the chip is transferred back to the MCU and steered through the USB port to a host computer where data are further analyzed with MATLAB. Neural recording signals from the database have to be preprocessed prior to being used as stimuli in the described platform. This involves the following steps: (i) bandpass-filtering for selecting the frequency band of interest (fourth-order Butterworth filters have been used), (ii) resampling for matching the frequency rates of the recordings and the ASIC (only if needed), (iii) adding noise for a given input-referred SNR value, (iv) digitizing the samples to the resolution of the ASIC (10-b), and (v) casting data in CSV format files for storing in the stick memory. As mentioned, files as large as 10 GB have been handled without any observed problem. B. Epilepsy Detection Results Differentblocksofintracranialdatafromthereferreddatabase were used for assessing the epilepsy detection capabilities of the DDA chip. These blocks were selected for their high confirmed seizure rate. No other considerations were taken into account. For each block, signals were first pre-filtered in MATLAB into conventional neuro-physiological bands. After a preliminary computer analysis of the records, the frequency band with the largest synchronization index activity was finally selected. For each experiment, the set of valid values were plotted alongside with the theoretical PLV value obtained through the Hilbert Transform method. For illustration purposes, Fig. 12 shows the synchronization results obtained from two channels of a recording block measured in a 48 years old, female patient. In this block, one low amplitude fast activity (lafa) seizure was annotated (dotted line). Before driving the ASIC, Gaussian noise were added to the signals for a spot SNR of 30 dB at 30 Hz. The plot shows the experimental synchronization index (labeled as “ASIC SI”) and the calculated PLV value (labeled as “SIM PLV”). The DDA algorithm detects all major changes in synchrony and, indeed, follows a similar trend as the PLV algorithm. This is verified by the high correlation measure of approximately 91% (in line with Fig. 7). Fig.12alsoshowsthattheseizuremanifestswitharemarkable increase in the synchronization measures. This feature is in fact on the basis of the detection mechanism presented in [20] and [28] in which seizures are detected by applying thresholds on the DELGADO-RESTITUTO et al.: PHASE SYNCHRONIZATION OPERATOR FOR ON-CHIP BRAIN FUNCTIONAL CONNECTIVITY COMPUTATION 965 TABLE V RECORDS USED IN THE VERIFICATION OF THE PHASE SYNCHRONIZATION PROCESSOR a(S/D/M/G) for surface, depth, strip and grid electrodes. Fig. 12. Experimental synchronization index and simulated PLV value from recordings at positions M5 and M8 of a 1×8EcoG strip implanted in patient FR_1084. Signals were pre-filtered in the β-band. The input rate of the chip was adjusted to the sample rate of the recordings, i.e., 1024 S/s. The Sdda curve (labeled “ASIC SI”) was obtained at a sample rate of 1 S/s. PLV indicators calculated between pairs of channels. In some cases, as in Fig. 12, seizures are preceded by a sudden drop in synchronization [18], however, we have not verified this feature in all the analyzed blocks. In this work, a similar thresholding approach has been followed and, as in the previous experimental results, PLV values were also calculated for comparison purposes. To this end, a simulation experiment with 10 recording sets included in the European Epilepsy database has been conducted. In total, 1803.2 hours of recordings and 61 annotated seizures have been analyzed. Simulations were carried out in MATLAB using mathematical models of the DDA and the PLV processors. In the DDA case, the model closely follows the block diagram of Fig. 4 while, in the PLV case, the model uses the architecture presented in [28]. In both cases, an input signal resolution of 10-b has been assumed. Table V shows the records used for the verification of the proposed phase synchronization calculation approach. For each patient the table shows the electrode configuration (both EEG and invasive), sampling rate, number of seizures, recording duration and dominant seizure pattern. The EEG data were acquired using a Neurofile NT digital video EEG system, while the invasive measurements were obtained with depth, strip or grid electrodes from Ad-Tech. For each patient, three different phase synchronization indexes, involving four intracranial electrodes, were calculated. In all cases, the three contact pairs have one electrode in common, located in the proximity of the epileptic focus, while the other three electrodes, also close to the position at which the seizure originates although not necessarily from the same array, were used as reference. These electrodes are identified in Table V as #F-{#R1, #R2, #R3}, where #F denotes the electrode at the seizure focus and #Rx denote the references. Prior to the calculation of the synchronization indexes, records were bandpass filtered in the frequency bands where higher activities wereobserved.Thisgavetheselectionof β-bandforall patientsbut patient #4,forwhichtheα-bandwasmostinsightful. In the case of patient #10, the β- and θbands offered similar results. In order to estimate the sensitivity (number of seizures which arecorrectlydetected dividedby the totalnumber of real positive cases) of the DDA and the PLV approaches in terms of the False Positive Rate, FPR(number of wrong positive detections along the interictal periods), different thresholds per contact pair and patient have been applied over the corresponding synchronization indexes. Thresholds for the DDA and the PLV algorithms do not necessarily coincide even if applied on the same electrode pair. In our experiment, a detection is produced whenever one out of three synchronization indexes (or PLV values) raise above theircorrespondingthresholds.Hence,decisionsare takenbased on the outcomes from three contact pairs, not just one. Additionally, the following considerations have been adopted [2]. An observation window of 15 min before the annotated seizure has been defined for the detection of true positives. Hence, if a crossing is detected any time within this window, which can be associatedtoapreictalstate,atrueseizuredetectionisaccounted for. Further, in order to exclude effects from postictal states, recording periods within 30 min after the onset of a seizure were discardedfromthresholdcalculations. If twosuccessive seizures are separated by less than 45 min, the preictal window before the onset of the following seizure is preserved in the analysis. Any period of time excluding the ictal events and the pre/post ictal paddings is regarded as interictal. No cool-off period after false positive alarms has been defined and, hence, all the wrong detections occurring in the interictal periods contribute to the FPR