Algorithm and software to automatically identify latency and amplitude features of local field potentials recorded in electrophysiological investigation
 Maria Rubega^{1},
 Claudia Cecchetto^{1, 2},
 Stefano Vassanelli^{2} and
 Giovanni Sparacino^{1}Email author
DOI: 10.1186/s1302901700625
© The Author(s). 2017
Received: 30 July 2015
Accepted: 28 January 2017
Published: 7 February 2017
Abstract
Background
Local field potentials (LFPs) evoked by sensory stimulation are particularly useful in electrophysiological research. For instance, spike timing and current transmembrane current flow estimated from LFPs recorded in the barrel cortex in rats and mice are exploited to investigate how the brain represents sensory stimuli. Recent improvements in microelectrodes technology enable neuroscientists to acquire a great amount of LFPs during the same experimental session, calling for algorithms for their quantitative automatic analysis. Several computer tools were proposed for LFP analysis, but many of them incorporate algorithms that are not open to inspection or modification/personalization. We present a MATLAB software to automatically detect some important LFP features (latency, amplitude, timederivative value in the inflectionpoint) for a quantitative analysis. The software features can be customized by the user according to his/her personal research needs. The incorporated algorithm is based on PhillipsTikhonov regularization to deal with noise amplification due to illconditioning. In particular, its accuracy in the estimation of the features of interest is assessed in a Monte Carlo simulation mimicking the acquisition of LFPs in different SNR (signaltonoiseratio) conditions. Then, the algorithm is tested by analyzing a real set of 2500 LFPs recorded in rat after whisker stimulation at different depths in the primary somatosensory (S1) cortex, i.e., the region involved in the cortical representation of touch in mammals.
Results
Automatic identification of LFP features by the presented software is easy and fast. As far as accuracy is concerned, error indices from simulated data suggest that the algorithm provides reliable estimates . Indeed, results obtained from LFPs recorded in rat after whisker stimulation are in line with the known sequential activation of the microcircuits of the S1 cortex.
Conclusion
A MATLAB software implementing an algorithm to automatically detect the main LFPs features was presented. Simulated and real case studies showed that the employed algorithm is accurate and robust against measurement noise. The available code can be used as it is, but the reported description of the algorithms allows users to easily modify the code to cope with specific requirements.
Keywords
LFP Automated analysis Rat barrel cortex Whisker stimulation PhillipsTikhonov regularization NeuroscienceBackground
Local field potentials (LFPs) reflect the synchronized population activity of several neurons recorded by smallsize electrodes in the brain. Any kind of transmembrane current in brain cells contributes to the extracellular fields known as LFPs. Thus, the amplitude and frequency of LFPs depend on the proportional contribution of the multiple sources and various properties of the brain tissue and its cells [1]. A quantitative analysis of LFPs is important in many neuroscience investigations. For instance, an interesting case study concerns how rats and other nocturnal animals process information about the spatial coordinates of objects and their identity by seeking out and palpating objects with their whiskers [2], a paradigm believed to be useful for studying how the brain represents sensory stimuli, also in human beings [3]. In particular, experiments in awake rats have demonstrated that the barrel cortex processes the whisking signal [4], while in anesthetized animals millisecondscale S1 firing patterns encode whisker stimuli [5]. In such studies, LFP signals are usually recorded from a barrel column of the rat S1 cortex using neural probes with multiple recording sites at different depths [6].
Visual identification and analysis of LFP features are obviously time consuming and exposed to the risk of subjectivity. Nowadays, a large number of LFPs (e.g. 500) can be collected within a single experimental session thanks to improvements in microelectrodes technology, making the development of algorithms for LFP quantitative automatic analysis even more urgent. Indeed, a number of algorithms and software tools for LFP analysis is available. For instance, in [8, 9], the software returns the location of the peaks at which data exceed an imposed threshold. In [10], spikes are first detected by setting an amplitude threshold, then are clustered by exploiting wavelet coefficients [12] or PCA (principal components analysis). In another software program [11], a graphical application for manual spike sorting is provided.
A limitation of many software tools for LFP analysis is that algorithmic details are not fully described and/or codes are not open to customization according to specific investigational needs. In the present paper, we illustrate a software to automatically detect some features of LFP waveforms pointed out in Fig. 1 named first maximum, inflection point and negative peak. The timederivatives used by the core algorithm to determine the features are calculated by the PhillipsTikhonov regularization method which mitigates measurement noise amplification due to illconditioning [13]. Several parameters of the MATLAB software can be customized by the user for personal research needs. The accuracy of the algorithm in correctly estimating the mentioned LFP features is first assessed in a Monte Carlo simulation mimicking the acquisition of LFP in different conditions of SNR. Then, the algorithm is tested by analyzing a set of real LFPs recorded in the S1 cortex of a rat at five different depths after whisker mechanical stimulation.
Methods
Brief description of the algorithm
With reference to the ideal case (noise free) of Fig. 1, the main goals are to detect: latency and amplitude of the first maximum (identified as the stimulus onset in this work), latency and amplitude of the negative peak, and first derivative value of the inflection point. In real world, recorded LFP are affected by measurement error and the automatic determination of their minima, maxima, and inflection points is complicated by the presence of noise which is amplified when timederivatives are computed because of illconditioning. This makes the straightforward use of techniques employed in other contexts difficult, e.g. for the automatic analysis of latencies of eventrelated potentials [14, 15].
 1.The singular value decomposition of matrix H = GF ^{−1} is performed to obtain Nsize unitary matrices U and V and matrix D = diag ([d_{1}, d_{2}, …, d_{N}]^{T})$$ {\mathbf{U}}^{\mathrm{T}}\mathbf{H}\mathbf{V}=\mathbf{D} $$(3)
This operation, which is performed only once, has numerical complexity O (N^{3}).
 2.γ is tuned until the condition of the discrepancy regularization criterion is met, i.e., until the following equation$$ \mathrm{WRSS}={\displaystyle {\sum}_{i=1}^N{\left(\frac{\gamma {\upxi}_{\mathrm{i}}}{{d_i}^2+\gamma}\right)}^2} $$(4)
is (approximately) matched, where ξ = U ^{T} y and WRSS stands for weighted residual sum of squares. As explained in [16], for each trial value of γ, the numerical complexity of this stage is only O (N).
 3.As γ is determined, the vector of the unknown timederivative samples û is computed as$$ \hat{\mathbf{u}}={\mathbf{F}}^{1}\mathbf{V}\mathbf{n} $$(5)
where n = [n_{1}, n_{2}, …, n_{i}, …, n_{N}] with \( {\mathrm{n}}_{\mathrm{i}}=\frac{{\mathrm{d}}_{\mathrm{i}}\;\upxi \mathrm{i}}{{{\mathrm{d}}_{\mathrm{i}}}^2+\upgamma} \). This stage has O (N^{3}) complexity.
As byproduct of the procedure, a regularized version of the LFP is obtained by multiplying û by matrix G. The difference between y and Gû over σ is the vector of the normalized residuals.
After the regularization step, the unknown first maximum and negative peak of the LFP are determined by finding where the estimated first timederivative crosses the zero line (change in the sign is also required). To make their estimation more robust avoiding the possibility to estimate a local minimum instead of the global one, a minimal distance (determined by the experimenter) between them is imposed. After estimating the first maximum and the negative peak, the onset is computed between them (the relative position of the onset is also forced by the user). Finally, the inflection point is estimated by finding where the second timederivative goes to zero in the interval between the previously determined first maximum and the negative peak (presence of change in its sign in the neighborhood interval is also verified).
MATLAB implementation

■ The path of the folder containing the data.

■ The name of the .mat file containing the data.

■ The name of the experiment (used to name the .mat file and .xls file that will contain the extracted features).

■ The depth of recording in order to name the .xls sheet (thus the .xls file will contain a sheet for each depth of recording referred to the same experiment).

■ The sampling frequency of the data.

■ The starting and the finishing time of the window of analysis.

■ The factor n to possibly decrease the sampling rate (if the sampling frequency of the acquired data is unnecessarily high for the purpose).

■ The position of the onset between the first maximum and the main negative peak, expressed in the range [0 1]. A latency parameter is computed by subtracting the stimulus onset to the negative peak time.

■ The minimum distance hypothesized by the experimenter from the onset and the main negative peak.
The .mat file containing the data has to be organized as follows: a matrix, in which each column stands for a single recording, and a time vector in ms (Additional file 1). If the data are not in the .mat format, a script to convert .txt file to .mat file is provided. The .txt file has to be structured in columns, in which the first one contains the time vector and the others the amplitudes of the recordings (Additional file 2; Additional file 3). Eventually, the software displays the results, and produces a .mat file and a .xls file containing the features extracted. This approach lets the user import the results into MATLAB or in Microsoft Excel for a further offline processing.
Results
Software use in a representative case study
To allow more flexibility than that allowed by the GUI, the software was intentionally organized in nested functions that investigators with basic knowledge of MATLAB can easily customize. Indeed, the main script (Additional file 4) calls two different functions (Additional file 5; Additional file 6), where the signal smoothing and the computation of its first and second timederivatives are performed following the stages described in the Methods section (see also Additional file 7; Additional file 8) . Thus, it is simple to do some little tweaks, e.g., changing the value of σ, and/or some bigger tweaks, e.g., modifying the criterion for the estimation of γ.
Reliability of the algorithm and robustness against SNR (simulated data)
Error indices in simulated data
Error indices  

SNR  t_{max} = t_{onset} [ms]  A_{max} = A_{onset} [unitless]  t_{peak} [ms]  A_{peak} [unitless]  1^{st} derivative value in the inflection point [unitless] 
10  0.25 (0.12)  −0.01 (0.14)  −0.16 (0.09)  0.01 (0.01)  0.05 (0.02) 
5  0.89 (0.96)  −0.01 (0.31)  −0.64 (0.36)  0.03 (0.02)  −0.21 (0.36) 
3  2.77 (1.24)  −0.73 (0.99)  −1.39 (1.09)  0.01 (0.03)  −0.06 (0.39) 
Overall results suggest that the method is sufficiently precise for a reliable automated detection of the LFP main features after whisker stimulation. Algorithm numerical efficiency is acceptable, e.g., execution time to process the 100 LFP recordings resulted 714 ms using an Intel Core i74790 at 3.6 GHz, with 16 GB of RAM.
Test on experimental data
Main features automatically detected in real data
Main features automatically extracted  

Depth [μm]  t_{max} = t_{onset} [ms]  A_{max} = A_{onset} [mV]  t_{peak} [ms]  A_{peak} [mV]  1^{st} derivative value in the inflection point [mV/ms] 
320 (layer II)  10.92 (0.09)  0.086 (0.002)  19.80 (0.11)  0.607 (0.010)  −0.038 (0.001) 
420 (layer III)  9.43 (0.04)  0.035 (0.002)  18.26 (0.06)  1.070 (0.015)  −0.064 (0.002) 
520 (layer IV)  8.76 (0.11)  0.022 (0.004)  18.75 (0.16)  1.005 (0.019)  −0.064 (0.001) 
720 (layer IV)  8.04 (0.04)  −0.008 (0.001)  17.29 (0.05)  1.122 (0.019)  −0.069 (0.002) 
920 (layer Va)  8.20 (0.30)  −0.006 (0.005)  18.2 (0.3)  0.764 (0.018)  −0.026 (0.001) 
Conclusions
In this paper, we presented a MATLAB software, exploiting PhillipsTikhonov regularization to automatically detect the first maximum, the following negative peak, and the first timederivative value in the inflection point between them, in LFPs evoked by whisker stimulation in rat barrel cortex. Preliminary experimental tests performed on simulated data with an increasing level of noise proved that the algorithm can be successfully used in automatically estimating the features of interest. Moreover, the method was successfully exploited to analyse a large LFP dataset to evaluate the differences in the response at different recordings depths. The codes were designed to allow customization by the user. Thus, on one hand, the investigator without programming background can customize the analysis by easily setting the main parameters to guide the features estimation and visualize its results, on the other hand, the programming user can modify the code exploiting its nested organization.
Current activity regards the massive analysis of large LFP datasets and their physiological interpretation. Moreover, the algorithm may be extended to estimate other features of interest, such as the duration of the positive deflection and of the negative valley of Fig. 1.
Notes
Abbreviations
 A_{max} :

Amplitude of the first maximum
 A_{onset} :

Amplitude of the onset
 A_{peak} :

Amplitude of the negative peak
 Inf:

Infinity
 LFPs:

Local field potentials
 P:

Postnatal day
 PCA:

Principal Components Analysis
 S1:

Primary somatosensory cortex
 sd^{2} :

Variance
 SNR:

Signaltonoise ratio
 t_{max} :

Latency of the first maximum
 t_{onset} :

Latency of the onset
 t_{peak} :

Latency of the negative peak
Declarations
Acknowledgements
Not applicable.
Funding
This study was supported by the University of Padova and from the 7th Framework Programme of the European Commission through “RAMP” project (http://www.rampproject.eu) with Contract No. 612058.
Availability of data and materials
Data and material supporting the results reported in the article can be found in http://bio.dei.unipd.it/rubega_scbm/.
Authors’ contributions
MR designed the algorithms, the software and the simulation; CC collected and analysed the data, and provided the results interpretation; SV supervised the design of the study; GS supervised the design of the algorithms. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
All experimental procedures were approved by the institutional animal welfare committee of the University of Padova, and were carried out in strict adherence to the Italian regulations on animal protection and care and with the explicit approval of the Italian animal welfare regulations (Decreto autorizzativo 447/2015PR).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Einevoll GT, Kayser C, Logothetis NK, Panzeri S. Modelling and analysis of local potentials for studying the function of cortical circuits. Neuroscience. 2013;14:770–85.PubMedGoogle Scholar
 Diamonds ME, von Heimendahl M, Knutsen PM, Kleinfeld D, Ahissar E. ‘Where’ and ‘what’ in the whisker sensorimotor system. Nat Rev Neurosci. 2008;9:601–12.View ArticleGoogle Scholar
 Knutsen PM, Pietr M, Ahissar E. Haptic object localization in the vibrissal. J Neurosci. 2006;26:8451–64.View ArticlePubMedGoogle Scholar
 Zuo Y, Safaai H, Notaro G, Mazzoni A, Panzeri S, Diamond ME. Complementary Contributions of Spike Timing and Spike Rate to Perceptual. Curr Biol. 2015;25:357–63.View ArticlePubMedGoogle Scholar
 Arabzadeh E, Panzeri S, Diamond ME. Deciphering the spike train of a sensory neuron: counts and temporal patterns in the rat whisker pathway. J Neurosci. 2006;26:9219–26.View ArticleGoogle Scholar
 Pettersen KH, Devor A, Ulbert I, Dale AM, Einevoll GT. Currentsource density estimation based on inversion of electrostatic forward solution: effects of finite extent of neuronal activity and conductivity discontinuities. J Neurosci Methods. 2006;154:116–33.View ArticlePubMedGoogle Scholar
 Kublic E. Contextual impact on sensory processing at the barrel cortex of awake rat. Acta Neurobiol Exp. 2004;64:229–38.Google Scholar
 Bokil HS, Andrews P, Kulkarni JE, Mehta S, Mitra P. Chronux: a new platform for analyzing neural signals. J Neurosci Methods. 2010;192:146–51.View ArticlePubMedPubMed CentralGoogle Scholar
 Bologna L, Pasquale V, Garofalo M, Baljon PL, Maccione A, Martinoia S, Chiappalone M. Investigating neuronal activity by SPYCODE multichannel data analyzer. Neural Netw. 2010;23:685–97.View ArticlePubMedGoogle Scholar
 Known KY, Eldawlatly S, Oweiss K. NeuroQuest: a comprehensive analysis tool for extracellular neural ensemble recordings. J Neurosci Methods. 2012;204:189–201.View ArticleGoogle Scholar
 Hazan L, Zugaro M, Buzsaki G. Klusters, NeuroScope, NDManager: a free software suite for neurophysiological data processing and visualization. J Nerosci Methods. 2006;155:207–16.View ArticleGoogle Scholar
 Quian Quiroga R, Nadasdy Z, BenShaul Y. Unsupervised spike detection and sorting with wavelets and superparamagnetic Clustering. Neural Comput. 2004;16:1661–87.View ArticlePubMedGoogle Scholar
 Tikhonov AN. Solution of incorrectly formuled problems and the regularization method. Soviet Math Dokl. 1963;4:1624.Google Scholar
 D’Avanzo C, Schiff S, Amodio P, Sparacino G. A Bayesian method to estimate singletrial eventrelated potentials with application to the study of the P300 variability. J Neurosci Methods. 2011;198:114–24.View ArticlePubMedGoogle Scholar
 Blankerts B, Lemm S, Treder M, Haufe S, Muller KR. Singletrial analysis and classification of ERP components – a tutorial. Neuroimage. 2011;56:814–25.View ArticleGoogle Scholar
 De Nicolao G, Sparacino G, Cobelli C. Nonparametric Input Estimation in Physiological Systems: Problems, Methods, and Case Studies. Automatica. 1997;33:851–70.View ArticleGoogle Scholar
 Twomey S. The application of numerical filtering to the solution of integral equations of the first and encountered in indirect sensing measurements. J Franklin Inst. 1965;279:95–109.View ArticleGoogle Scholar
 Frigo, G., Rubega, M., Lezziero, G., Fontana, R., Cecchetto, C., Vassanelli, S., Sparacino, G. & Bertocco, M. (2015, May). A softwarebased platform for multichannel electrophysiological data acquisition. In Medical Measurements and Applications (MeMeA), 2015 IEEE International Symposium on (pp. 353358). IEEE.
 Fox K. Barrel Cortex. Cambridge: Cambridge University Press; 2008.View ArticleGoogle Scholar
 Jellema T, Brunia CHM, Wadman WJ. Sequential activation of microcircuits underlying somatosensoryevoked potentials in rat neocortex. Neuroscience. 2004;129:283–95.View ArticlePubMedGoogle Scholar
 Mahmud M, Pasqualotto E, Bertoldo A, Girardi S, Maschietto M, Vassanelli S. An automated method for detection of layer activation order in information processing pathway of rat barrel cortex under mechanical whisker stimulation. J Neurosci Methods. 2011;196:141–50.View ArticlePubMedGoogle Scholar
 Swanson LW. Brain Maps: Structure of the Rat Brain. London: Academic; 2003.Google Scholar