Neural Networks 68 (2015) 89–95
Contents lists available at ScienceDirect
Neural Networks journal homepage: www.elsevier.com/locate/neunet
The transfer function of neuron spike Igor Palmieri a,∗ , Luiz H.A. Monteiro b,a , Maria D. Miranda a a
Departamento de Engenharia de Telecomunicações e Controle, Escola Politécnica da Universidade de São Paulo, Av. Prof. Luciano Gualberto, travessa 3, n. 158, 05508-900, São Paulo, SP, Brazil b Escola de Engenharia da Universidade Presbiteriana Mackenzie, Pós-graduação em Engenharia Elétrica, Rua da Consolação, n. 896, 01302-907, São Paulo, SP, Brazil
article
info
Article history: Received 20 December 2014 Received in revised form 9 March 2015 Accepted 8 April 2015 Available online 22 April 2015 Keywords: Action potential Neuronal spike Parametric modeling Transfer function
abstract The mathematical modeling of neuronal signals is a relevant problem in neuroscience. The complexity of the neuron behavior, however, makes this problem a particularly difficult task. Here, we propose a discrete-time linear time-invariant (LTI) model with a rational function in order to represent the neuronal spike detected by an electrode located in the surroundings of the nerve cell. The model is presented as a cascade association of two subsystems: one that generates an action potential from an input stimulus, and one that represents the medium between the cell and the electrode. The suggested approach employs system identification and signal processing concepts, and is dissociated from any considerations about the biophysical processes of the neuronal cell, providing a low-complexity alternative to model the neuronal spike. The model is validated by using in vivo experimental readings of intracellular and extracellular signals. A computational simulation of the model is presented in order to assess its proximity to the neuronal signal and to observe the variability of the estimated parameters. The implications of the results are discussed in the context of spike sorting. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction The mathematical modeling of neuronal signals is one of the most relevant problems in neuroscience. This question has been driven by technological advances that allowed the experimental acquisition of neuronal signal in several situations, establishing the basis for many quantitative models found in literature (Koch & Segev, 1998; Sterratt, Graham, Gillies, & Willshaw, 2011). Such models have been widely employed to solve computational tasks in a computer (Gao, Yang, Cai, & Liu, 2012; Menezes & Monteiro, 2011) and to understand the functioning of nervous systems (Monteiro, Bussab, & Chaui-Berlinck, 2002; Steuber & Jaeger, 2013). One of the first models was proposed by Huxley and Hodgkin in the 1950s, using a system of continuous-time equations based on experiments with giant squid cells (Hodgkin & Huxley, 1952). This model served as inspiration for subsequent proposals, mainly focused on the reproduction of the action potential dynamics, as in FitzHugh (1961). Discrete-time versions (Ibarz, Casado, & Sanjuán, 2011), as the McCulloch–Pitts model, can be considered simplifications to facilitate simulations in digital environments.
∗
Corresponding author. E-mail addresses: [email protected] (I. Palmieri), [email protected] (L.H.A. Monteiro), [email protected] (M.D. Miranda). http://dx.doi.org/10.1016/j.neunet.2015.04.003 0893-6080/© 2015 Elsevier Ltd. All rights reserved.
There are many biophysical and chemical factors involved in the generation of neuronal signals; hence, to model such signals is a particularly difficult task. Moreover, the neuron behavior is highly influenced by the connections with surrounding nerve cells, forming a network which itself poses a great challenge for mathematical representation (Dayan & Abbott, 2005). Most of the referred models are designed with the main intent of representing the neuronal dynamics and the interaction with external agents. However, in several situations, a parametrization associated solely with the signal waveform is of great interest. Extra and intracellular action potential recordings have been used to tune the parameter values of continuous-time nonlinear models, which explicitly consider the physical properties of the neuronal membrane. This tuning is commonly performed by using parameter estimation techniques (Gold, Henze, Koch, & Buzsáki, 2006; Tabak, Murphey, & Moore, 2000; Toth, Kostuk, Meliza, Margoliash, & Abarbanel, 2011; Wang et al., 2014; Willms, Baro, Harris-Warrick, & Guckenheimer, 1999). However, the number of parameters to be matched can be as great as 100 (Gold, Henze, & Koch, 2007). In fact, the complexity of these models can impair computational simulations and analytical studies involving many neurons. There are also discrete-time nonlinear models that can generate a signal which is qualitatively similar to neuronal recordings (Aihara, Takabe, & Toyoda, 1990; GirardiSchappo, Tragtenberg, & Kinouchi, 2013; Hsu, Gobovic, Zaghloul,
90
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
& Szu, 1996; Mesbah, Moghtadaei, Golpayegani, & Towhidkhah, 2014). However, these models neither take into account details of the neurophysiological processes nor experimental data. Their parameter values are empirically set. The parametrization of intracranial signals is of great interest in the context of the identification problem. The study of efficient techniques to isolate the firings of each active cell from an extracellular mixture containing sequences of spikes from several neurons is a problem extensively discussed in Neuroscience, formally known as spike sorting. Usually these techniques are implemented in three stages (Gibson, Judy, & Markovic, 2012; Lewicki, 1998; Quiroga, Nádasdy, & Ben-Shaul, 2004; Rutishauser, Schuman, & Mamelak, 2006): from raw electrical potentials recorded by using arrays of electrodes, spikes are detected, then parameterized, and finally sorted, attributing every single spike observed to a particular neuron. In this context, a relevant question is how the parameter sets describe effectively the neuronal activities over time. A simplified model of how the signal that arrives at the electrodes is formed from a stimulus and its evolution over time can possibly help to elucidate this question. The main goal of this paper is to propose a discrete-time linear time-invariant (LTI) model with a rational function, to represent the neuronal spike observed by an electrode located in the surroundings of the nerve cell. This approach provides a low-complexity alternative to the neuronal spike model, employing system identification and signal processing concepts. It is dissociated from any considerations about the biophysical processes of the neuronal cell, and it is not directly derived from any traditional neuronal models. The model is validated by using in vivo experimental readings of intracellular and extracellular signals. Although the proposed model does not fully comprise the cellular dynamics, it allows a straightforward simulation of the neuronal activity by a pair of linear equations, once its parameters are obtained. The remainder of this paper is organized into five sections. In Section 2, the experimental methods involved in the data acquisition of neuronal electrical signals are briefly described. The proposed model is introduced in Section 3, in the form of a cascade association of two subsystems: one that generates an action potential from an input stimulus, and one that represents the medium between the cell and the electrode. The parameters for these subsystems are estimated from experimental data in Section 4. In Section 5, a computational simulation of the model is presented in order to assess its proximity to the neuronal signal and the variability of the estimated parameters. The implications of the results are discussed in Section 6, considering the limitations of the proposed approach. 2. Experimental methods The estimation of the proposed model parameters is performed by using in vivo experimental data recordings of neuronal signals generated by mice, provided by CRCNS (Collaborative Research in Computational Neuroscience) (Henze, Harris et al., 2000). For each experiment in the database, there are time series of action potentials obtained directly from the cell body through glass micropipettes, and of the electric potential observed simultaneously in the external environment, through electrodes placed in an area near the chosen cell. These external electrodes are mounted in arrays of four to six uniformly spaced elements, providing spatial diversity in the recordings of the electric potential outside the cell. All the recordings were obtained from the region of the hippocampus; a great part of the pulses corresponds to the pyramidal cell activity observed in a layer known as CA1. The time series were digitized at 12 or 16 bits, and sampled at frequency
Fig. 1. (a) Segment of the original acquired extracellular signal; (b) the same segment after filtering to remove low frequency components.
of 10, 20, 25 or 50 kHz, depending on the experiment. The values were then recorded in binary files, and made publicly available by CRCNS at their website (Henze, Harris et al., 2000). According to the coordinators of the experiment, the extracellular signal was first amplified 1000–8000 times, and filtered at the source by a low pass filter with cutoff frequency of 3 kHz (Henze, Borhegyi et al., 2000). Before using the digitalized data in the modeling process, however, we also employed a filtering process in each sequence. In the extracellular sequence, we applied an elliptical-type high-pass filter with cutoff frequency of 300 Hz. The purpose of this filtering is to remove low frequency components, mainly associated with the LFP (Local Field Potential). Although it also contains some information about the neuronal activity (Scherberger, Jarvis, & Andersen, 2005), this fact is not addressed in the proposed parametrization. Fig. 1 shows a fragment of the extracellular signal before and after this stage. The intracellular sequence, composed almost exclusively by the action potential of the neuron, is also passed through a high-pass filter, but with cutoff frequency of 50 Hz. This filtering only aims to remove low frequency components, especially in zero frequency, to ensure the level of the resting potential near zero. After the filtering process, we selected a fixed interval of 25 s within one of the available recordings, containing 86 action potentials of the same neuron. We noted that the corresponding intracellular series contains solely these action potentials, as it is obtained directly in the membrane of a single neuron, while the extracellular series contains the spikes of all neurons in the region during the same time interval. Centered around the exact sample of the maximum value of each signal, we defined a window of 14 ms (N = 140 samples for a sampling rate of 10 kHz). In the remainder of this paper, the intracellular and extracellular signals inside each time window are denoted as {s(ℓ)} and {y(ℓ)}, respectively, and the estimate of these signals are denoted as {ˆs(ℓ)} and {ˆy(ℓ)} with ℓ = 0, . . . , N − 1. 3. The proposed model The main goal of parametric modeling is to build a mathematical model able to represent relevant characteristics of a given signal by a reduced set of parameters. Here we consider the specific case of modeling the discrete-time extracellular signal y(n) whose sample set is obtained through direct observation inside a fixed time window. Fig. 2 shows the simplified block diagram of the proposed model, in which x(n) is the input and the output yˆ (n) is the estimate of y(n). We assume that this signal can be represented by a LTI
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
91
Fig. 2. Simplified block diagram of the proposed model: F (z ) is the all-pole model, G(z ) is the moving average model, yˆ (n) is the output signal, and η(n) is an additive noise include to yield a signal similar to the recordings.
model with transfer function H (z ). This function is represented with two subsystems F (z ) and G(z ) associated in series, so that q
gm
b k z −k
k=0 p
H (z ) , F (z ) G(z ) =
1−
(1) ak z − k
k=1
in which gm is a gain factor, and p and q represent the orders of the polynomials at the denominator and at the numerator, respectively. The parametrization process corresponds to estimate H (z ), or the values of the coefficients ak and bk (and also the gain gm ), that together define a model capable of yielding an approximation yˆ (n) of the original signal y(n), given an input x(n). Essentially, when a punctual stimulus is above a threshold, the neuron fires. If the stimulus received by the neuron is below such a threshold, it can be considered null; if it is above, its value is taken as equal to one. This assumption is common in artificial neural networks, where the neuron output is usually a binary variable; typically, 0 or 1 (Haykin, 2009). As a premise, in our approach, it is considered that y(n) is basically a deterministic signal, and the input of the model is a Kronecker delta, or discrete unit impulse (Hayes, 1996), time shifted to the exact moment when the neuron fired. This unit impulse x(n) = δ(n − ∆) =
1, 0,
n=0 n ̸= 0
(2)
is essentially a punctual stimulus on the nerve cell enough to make it fire at n = ∆. So, in our approach, the estimate of the observed extracellular spike output yˆ (n) has the same role as the impulse response of the system. Hence, H (z ) can be considered the transfer function of this system. The subsystem F (z ) represents the model associated with the action potential produced in the cell membrane, and thus it represents characteristics inherent to each neuron. We assume that the input–output relation is an autoregressive all-pole model, containing the denominator of the transfer function, scaled by a gain factor gm , that is, gm
F (z ) = 1−
p
.
(3)
ak z − k
k=1
Then the output sˆ(n) generated by this model is calculated recursively by the relation sˆ(n) =
p k=1
ak sˆ(n − k) + gm x(n) = sˆ Tp (n − 1) ap + gm x(n)
(4)
in which ap = [a1 a2 . . . ap ]T and sˆ p (n − 1) = [ˆs(n − 1) sˆ(n − 2) . . . sˆ(n − p)]T are the coefficients and the input column vectors respectively, and T denotes the vector transpose operator. In broad terms, the output sˆ(n) should be a good estimate of the original intracellular signal s(n). If the signal s(n) could be entirely represented by an all-pole system, so a proper choice of parameters would result in sˆ(n) = s(n). The subsystem G(z ) represents the medium between the cell and the electrode, and it acts as a communication channel through which the generated signal propagates. We assume that the input–output relation is a moving average model, that is, G(z ) =
q
b k z −k .
(5)
k=0
Then, its output is calculated by the relation yˆ (n) =
q
bk sˆ(n − k) = sˆ Tq+1 (n) bq+1
(6)
k=0
in which bq+1 = [b0 b1 . . . bq ]T and sˆ q+1 (n) = [ˆs(n) sˆ(n − 1) . . . sˆ(n − q)]T are the coefficients and the input column vectors, respectively. It is desired that the output of the complete system, formed by the concatenation of F (z ) with G(z ), provides an estimate of the signal y(n), recorded from the outside medium of the cell, with characteristics that depend on the neuron itself and the environment shared by the cells in that specific brain region. The output of the system formed by the cascade of the autoregressive and moving average models comprises the deterministic part of the signal that we are trying to model. However, observations of actual extracellular signal indicate the presence of some level of noise η(n) due to the measurement process. This effect may be represented by adding a random signal at the output of the function of the model. This noise can also take in account the variability in the underlying neural circuits. 4. Parameter estimation of the transfer function There are several approaches that can be used to calculate the coefficients of the model in order to approximate the actual signal in each block of the unknown system (Ljung, 1999). This work considers Prony’s method (Hayes, 1996) to produce the estimates of the autoregressive coefficients, associated with the model F (z ) of intracellular signal. In the case of the model G(z ), the optimal coefficients are estimated by using a system identification configuration whose input is the intracellular signal s(n), and the output is the extracellular signal y(n). In both cases we apply the Least Squares Error (LSE) criterion. For each modeled action
92
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
Fig. 3. Average intracellular signal and limits defined by the standard deviation.
Fig. 4. Average extracellular signal and limits defined by the standard deviation.
potential and spike, only a finite number N of samples are available. In our approach for estimated all the coefficients of F (z ) and G(z ), we consider the autocorrelation method which makes the assumptions that the input data prior to n = 0 and the data after n = N − 1 are zero (Oppenheim & Schafer, 2009).
this in mind, we define the estimation error between the desired signal y(n) and model output yˆ (n) as ey (n) = y(n) − yˆ (n)
(10)
in which yˆ (n) = sˆ Tq+1 (n) bq+1 . The solution for the best transversal filter, according to LSE criterion, is given by the equation
4.1. The parameters of F (z )
b q +1 =
N +q
−1 sq+1 (ℓ)
sTq+1
(ℓ)
N +q
sq+1 (ℓ)y(ℓ).
(11)
Prony’s method is a way of estimating the coefficients of an autoregressive and moving-average model, especially useful when applied to the case of the deterministic all-pole model (Hayes, 1996). The coefficients are obtained by solving the problem of forward prediction of the desired signal. The forward prediction error is defined as
Here, the inverse of the autocorrelation matrix is estimated by using the experimental readings of the intracellular signal {s(ℓ)}, as in Eq. (9), but the cross-correlation vector is estimated by using the intracellular {s(ℓ)} and extracellular {y(ℓ)} signals.
es (n) = s(n) − sˆ(n)
5. Results and simulation
(7)
in which sˆ(n) = (n − 1)ap . The optimal coefficients are obtained in order to minimize the energy of es (n). A compact way to represent this criterion is sTp
ap = arg min ap
N −1
|es (n)|
2
(8)
n =0
where ap is the vector of optimal coefficients. The solution that satisfies (8), according to the LSE criterion, must also satisfy the relation
ap =
−1
N +p−1
sp (ℓ) (ℓ) sTp
N +p−1
(9)
ℓ=0
ℓ=0
sp (ℓ)s(ℓ + 1) .
(I )
(II )
It is worth to note that the inverse of the autocorrelation matrix (I ) and the crosscorrelation vector (II ) are estimated by using in vivo experimental readings of intracellular signals {s(ℓ)}. The value of gm is calculated in a later step, as a gain factor to equalize the energy of the observed signal and the estimated signal. The last task of intracellular signal modeling is to define exactly the input x(n). For this end, an integer delay value ∆m is determined empirically, to be applied on the input signal, such that x(n) = δ(n − ∆m ) acts as a trigger in the same instant that the real neuron fired. 4.2. The parameters of G(z ) Given the input s(n) and the output y(n) of an unknown system, the system identification configuration allows us to obtain the coefficients bk of the moving average model (Haykin, 1996). With
ℓ=0
ℓ=0
In the following, we present the selected firings contained in the filtered signals from the database. Then, we show how the orders of the subsystems F (z ) and G(z ) are estimated. Finally, all the parameters values are calculated. With these coefficients, we obtain the estimates sˆ(n) and yˆ (n), which are compared with the original signals s(n) and y(n). 5.1. In vivo experimental readings of intracellular and extracellular signals We selected a 25 s segment of the intra and extracellular signals filtered from the original time series of the database, containing a sequence of 86 action potentials. Each action potential is represented by a window with N = 140 samples (14 ms). Despite the fact that the firing sequence has been generated by the same neuron, the waveform of the action potential and the spike do not have a fixed shape over time. Figs. 3 and 4 illustrate the mean values of intracellular and extracellular signals inside the 14 ms window, considering all the firing instants in the chosen segment. Observe that there is a noticeable variability in the waveform of both signals, and also an apparent level of noise in the extracellular sequence. 5.2. The estimate of the model order For the subsystem F (z ) and G(z ), it is necessary to define a number of coefficients that allow the model to reproduce the signal of interest with reasonable assertiveness. In other words, we must find the integers p and q that suffice to build a suitable model. This task is directly linked with the complexity of the optimal coefficients estimations in Eqs. (9) and (11).
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
93
Fig. 5. Partial autocorrelation function of the intracellular signal, with a confidence interval based on the standard error (SE) to help the definition of the model order. Fig. 7. Model output sˆ(n) vs. acquired intracellular signal s(n).
Fig. 6. Model order vs. normalized (by the energy of the real signal) mean squared error. Fig. 8. Model output yˆ (n) vs. acquired extracellular signal y(n).
In the all-pole model case, a good indication of the number of coefficients of the F (z ) is given by the partial autocorrelation function of the signal being modeled, while its value is far from zero (Ljung, 1999). This measure, however, does not make sense to the moving average extracellular model G(z ). As an alternative, in both subsystems, it is also possible to calculate the normalized mean squared error between the model output and the observed signal as function of the model order (Oppenheim & Schafer, 2009). When this error approaches a constant level, there is an indication that the addition of coefficients is no longer relevant. Fig. 5 shows the average partial autocorrelation function calculated with the intracellular time series, and the 95% confidence interval of the standard error (Box, Jenkins, & Reinsel, 2008). It is also possible to simulate both models with a different number of coefficients and observe their residual error. Fig. 6 exhibits the normalized mean squared error of intracellular and extracellular signals (given by Eqs. (7) and (10), respectively) of a single neuron, compared to the model order. These curves together suggest that we can use p = 6 and q = 3 coefficients for the intra and extracellular models, respectively. 5.3. The estimate of the intra and extracellular signals The samples of each window are then used to estimate the autocorrelation and the crosscorrelation functions, which will solve Eqs. (9) and (11). Each window therefore corresponds to a single set of model parameters, for both the intracellular and the extracellular models. We can better observe the result of the modeling process by analyzing the time interval of a single isolated action potential and its corresponding spike. Fig. 7 illustrates the output produced by the all-pole intracellular model along with the signal used to estimate the coefficients. We can see that there are some differences in the generated signal, especially in the repolarization phase.
The moving average extracellular subsystem is also simulated within the same time frame. In Fig. 8, the output of the complete model is shown, compared to the observed extracellular spike. Note that there is an influence of the experimental noise in this signal, which accounts for the main differences from the model output. 5.4. Distributions of the coefficients The variability of the waveform in the intra and extracellular signals is reflected in the changes of the estimates of the F (z ) and G(z ) coefficients, respectively. Figs. 9 and 10 present the normalized histogram (such that the sum of the bars is equal to one) of the coefficients calculated for each subsystem using the chosen sequence of firings. It is possible to observe the behavior of the model coefficients for multiple firings in the same sequence by looking at the statistical distribution of their values. Another important result is that the average coefficients of the extracellular model, seen in Fig. 10, indicate a transfer function of a differentiator, as expected given the capacitive characteristics of the external medium through which the signal propagates (Gibson et al., 2012). By inserting these average coefficients values in Eq. (4), with the average gain, we obtain the following relation to calculate the output for the observed neuron: sˆ(n) = sˆ T6 (n − 1) a6 + 14.622 δ(n − ∆)
(12)
in which a6 = [2.436 − 2.430 1.456 − 0.714 0.267 − 0.045]T . These coefficients form a system whose poles are inside the unitary circle, and consequently this is a bounded-input, boundedoutput (BIBO) stable system. Considering that the neuron fired at the instant n = ∆, the resulting values for the output will start at zero and then follow the sequence sˆ(∆ − 1) = 0; sˆ(∆) = 14.622; sˆ(∆ + 1) = 35.619; sˆ(∆ + 2) = 51.237, and so on, decreasing asymptotically to the initial signal level.
94
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
types of neurons (Bean, 2007). Thus, although Eqs. (12) and (13) were obtained from a particular database, the method here proposed can be used for modeling spikes with distinct waveforms. This is relevant because neural communication is based not only on frequency and phase of action potentials, but also in their amplitude (Bean, 2007; de Polavieja et al., 2005; Zheng et al., 2009). 6. Discussion
Fig. 9. Distribution of the estimated coefficients of the intracellular model, with average values indicated.
The model proposed in this paper has a high level of abstraction regarding any physiological process and presents a considerable level of simplification, being comparable to most of the discretetime models found in the literature (Dayan & Abbott, 2005). Its main feature is to employ a linear equation for representing the neuronal signal during the firing of the cell, instead of a set of non-linear equations. The output signal of the model is readily calculated once its parameters are estimated. Usually, neuronal models try to reproduce or mimic the dynamics involved in the generation of an action potential, including features like fire rate or the cell’s membrane conductances (Ibarz et al., 2011; Izhikevich, 2007). In the proposed model, the goal is to describe mainly its waveform. The generation of the action potential and the resulting spike are triggered by the input signal, containing a unit impulse time shifted to the exact moment of the firing. Despite showing minor deviations from the original signal, the main significance of this work is the structure for the neuronal signal model: autoregressive in the intracellular and moving average in the extracellular. It should also be noted that the parameters of the transfer function H (z ) vary between each firing, reflecting the changes in the signal waveform. These variations are observed in both intra and extracellular signals, even though they present different signal-to-noise ratio. Expected developments
Fig. 10. Distribution of the estimated coefficients of the extracellular model, with average values indicated.
An analogous relation can be calculated for the extracellular model, using the average coefficients in the Eq. (6): yˆ (n) = sˆ T4 (n) b4
(13)
in which b4 = [−0.454 0.452 − 0.275 0.228] . This signal also departs from zero, and the input remains zero till the neuron fires. The output sequence will follow the values yˆ (∆ − 1) = 0; yˆ (∆) = −6.638; yˆ (∆ + 1) = −9.562; yˆ (∆ + 2) = −11.183, also decreasing back to the same initial values. The dispersion around their mean values is a consequence of the signal variability, as displayed in Figs. 3 and 4. Such changes in the shape of action potentials can be due to differences in ionic conductances (Djouhri, Bleazard, & Lawson, 1998); changes on the levels of neurohormones, neuromodulators, neurotransmitters (Rodríguez-Molina et al., 2014); diversity on the stimulus history (de Polavieja, Harsch, Kleppe, Robinson, & Juusola, 2005); distinct T
In practical applications, an electrode employed in neural signal acquisition registers the combined activity of a set of neighbor neurons. The goal of spike sorting processes is to distinguish the neurons generating the activity detected by every electrode. This task includes stages of clustering and classification of the detected spikes, and is highly based on a parametric representation of the extracellular signal (Gibson et al., 2012; Lewicki, 1998; Quiroga et al., 2004; Rutishauser et al., 2006; Sarinho Filho, Eisencraft, Suyama, Fonoff, & Miranda, 2011). The results obtained with the proposed approach suggest further investigation of the connection between the polynomial H (z ) and the parametrization methods of spike sorting. However, simulations showed that the changes in the spike waveform are directly reflected in the model coefficients, even when the presence of noise is ignored. Consequently, it should be expected that any method of spike identification based on the extracellular signals from experimental data will also subject to the effects of variations in waveform parameters. Due to changes in the waveform of the spike, errors in parametrization and classification phases of any spike sorting process are hard to overcome. This suggests that the dynamics of the neuronal activity can be better described if the decision about the neuron activation is observed over time consecutively with the detection stage and not in a posterior stage like is commonly done in spike sorting. In such context, the use of adaptive techniques combined with the detection and the parameter estimation stages is of high interest, as addressed in Miranda and Gerken (1997) and Miranda, Silva, and Nascimento (2008). A deeper investigation should explore the role of the model parameters in the process of neuron identification.
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
The model can also be further expanded by considering simultaneous readings of different extracellular electrodes, resulting in a system containing multiple output channels. Since usually the electrodes are located in different points of the external medium, the resulting model should reflect the spatial diversity of the extracellular signal. A more general structure should also include more than one input, each one associated with the firings of an active cell in a region containing several neurons. However, this model would require observations of the intracellular signal of other cells. Due to the availability of the experimental data, consisting of only one intracellular time series, the simulations in this work were so far restricted to the single input case. Acknowledgment L.H.A. Monteiro is partially supported by National Council for Scientific and Technological Development (CNPq) under the grant #305827/2014-6. References Aihara, K., Takabe, T., & Toyoda, M. (1990). Chaotic neural networks. Physics Letters A, 144, 333–340. Bean, B. P. (2007). The action potential in mammalian central neurons. Nature Reviews Neuroscience, 8, 451–465. Box, G. E. P., Jenkins, G. M., & Reinsel, G. C. (2008). Time series analysis: forecasting and control (4th ed.). Hoboken, NJ: John Wiley & Sons. Dayan, P., & Abbott, L. F. (2005). Theoretical neuroscience: computational and mathematical modeling of neural systems. Cambridge, MA: MIT Press. de Polavieja, G. G., Harsch, A., Kleppe, I., Robinson, H. P. C., & Juusola, M. (2005). Stimulus history reliably shapes action potential waveforms of cortical neurons. The Journal of Neuroscience, 25, 5657–5665. Djouhri, L., Bleazard, L., & Lawson, S. N. (1998). Association of somatic action potential shape with sensory receptive properties in guinea-pig dorsal root ganglion neurones. The Journal of Physiology, 513, 857–872. FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1, 445–466. Gao, D., Yang, Z., Cai, C., & Liu, F. (2012). Performance evaluation of multilayer perceptrons for discriminating and quantifying multiple kinds of odors with an electronic nose. Neural Networks, 33, 204–215. Gibson, S., Judy, J. W., & Markovic, D. (2012). Spike sorting: the first step in decoding the brain. IEEE Signal Processing Magazine, 29, 124–143. Girardi-Schappo, M., Tragtenberg, M. H. R., & Kinouchi, O. (2013). A brief history of excitable map-based neurons and neural networks. Journal of Neuroscience Methods, 220, 116–130. Gold, C., Henze, D. A., & Koch, C. (2007). Using extracellular action potential recordings to constrain compartmental models. Journal of Computational Neuroscience, 23, 39–58. Gold, C., Henze, D. A., Koch, C., & Buzsáki, G. (2006). On the origin of the extracellular action potential waveform: a modeling study. Journal of Neurophysiology, 95, 3113–3128. Hayes, M. H. (1996). Statistical digital signal processing and modeling. New York, NY: John Wiley & Sons. Haykin, S. (1996). Adaptive filter theory (3rd ed.). Upper Saddle River, NJ: Prentice Hall. Haykin, S. (2009). Neural networks: a comprehensive foundation. Upper Saddle River, NJ: Prentice Hall. Henze, D. A., Borhegyi, Z., Csicsvari, J., Mamiya, A., Harris, K. D., & Buzsáki, G. (2000). Intracellular features predicted by extracellular recordings in the hippocampus in vivo. Journal of Neurophysiology, 84, 390–400. Henze, D.A., Harris, K.D., Borhegyi, Z., Csicsvari, J., Mamiya, A., & Hirase, H. et al. 2000. Simultaneous intracellular and extracellular recordings from hippocampus region CA1 of anesthetized rats. CRCNS.org. URL: http://dx.doi.org/10.6080/K02Z13FP.
95
Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117, 500–544. Hsu, C. C., Gobovic, D., Zaghloul, M. E., & Szu, H. H. (1996). Chaotic neuron models and their VLSI circuit implementations. IEEE Transactions on Neural Networks, 7, 1339–1350. Ibarz, B., Casado, J. M., & Sanjuán, M. A. F. (2011). Map-based models in neuronal dynamics. Physics Reports, 501, 1–74. Izhikevich, E. M. (2007). Dynamical systems in neuroscience. Cambridge, MA: MIT Press. Koch, C., & Segev, I. (Eds.) (1998). Methods in neuronal modeling: from ions to networks (2nd ed.). Cambridge, MA: MIT Press. Lewicki, M. S. (1998). A review of methods for spike sorting: the detection and classification of neural action potentials. Network: Computation in Neural Systems, 9, 53–78. Ljung, L. (1999). System identification: theory for the user (2nd ed.). Upper Saddle River, NJ: Prentice Hall. Menezes, R. A., & Monteiro, L. H. A. (2011). Synaptic compensation on Hopfield network: implications for memory rehabilitation. Neural Computing and Applications, 20, 753–757. Mesbah, S., Moghtadaei, M., Golpayegani, M. R. H., & Towhidkhah, F. (2014). Onedimensional map-based neuron model: a logistic modification. Chaos, Solitons and Fractals, 65, 20–29. Miranda, M. D., & Gerken, M. (1997). A hybrid least squares QR-lattice algorithm using a priori errors. IEEE Transactions on Signal Processing, 45, 2900–2911. Miranda, M. D., Silva, M. T. M., & Nascimento, V. H. (2008). Avoiding divergence in the Shalvi–Weinstein algorithm. IEEE Transactions on Signal Processing, 56, 5403–5413. Monteiro, L. H. A., Bussab, M. A., & Chaui-Berlinck, J. G. (2002). Analytical results on a Wilson–Cowan neuronal network modified model. Journal of Theoretical Biology, 219, 83–91. Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-time signal processing (3rd ed.). Upper Saddle River, NJ: Prentice Hall. Quiroga, R. Q., Nádasdy, Z., & Ben-Shaul, Y. (2004). Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering. Neural Computation, 16, 1661–1687. Rodríguez-Molina, V., Patiño, J., Vargas, Y., Sánchez-Jaramillo, E., Joseph-Bravo, P., & Charli, J. L. (2014). TRH regulates action potential shape in cerebral cortex pyramidal neurons. Brain Research, 1571, 1–11. Rutishauser, U., Schuman, E. M., & Mamelak, A. N. (2006). Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo. Journal of Neuroscience Methods, 154, 204–224. Sarinho Filho, J. N. S., Eisencraft, M., Suyama, R., Fonoff, E. T., & Miranda, M. D. (2011). The use of least squares lattice algorithm in the parameterization and sorting of action potentials signals. In Bioelectronics, biomedical, and bioinspired systems V and nanotechnology V . 80680Q. Scherberger, H., Jarvis, M. R., & Andersen, R. A. (2005). Cortical local field potential encodes movement intentions in the posterior parietal cortex. Neuron, 46, 347–354. Sterratt, D., Graham, B., Gillies, A., & Willshaw, D. (2011). Principles of computational modelling in neuroscience. Cambridge, UK: Cambridge University Press. Steuber, V., & Jaeger, D. (2013). Modeling the generation of output by the cerebellar nuclei. Neural Networks, 47, 112–119. Tabak, J., Murphey, C. R., & Moore, L. (2000). Parameter estimation methods for single neuron models. Journal of Computational Neuroscience, 9, 215–236. Toth, B. A., Kostuk, M., Meliza, C. D., Margoliash, D., & Abarbanel, H. D. I. (2011). Dynamical estimation of neuron and network properties I: variational methods. Biological Cybernetics, 105, 217–237. Wang, R., Wang, J., Deng, B., Liu, C., Wei, X., Tsang, K. M., et al. (2014). A combined method to estimate parameters of the thalamocortical model from a heavily noise-corrupted time series of action potential. Chaos, 24, 013128. Willms, A. R., Baro, D. J., Harris-Warrick, R. M., & Guckenheimer, J. (1999). An improved parameter estimation method for Hodgkin–Huxley models. Journal of Computational Neuroscience, 6, 145–168. Zheng, L., Nikolaev, A., Wardill, T. J., O’Kane, C. J., de Polavieja, G. G., & Juusola, M. (2009). Network adaptation improves temporal representation of naturalistic stimuli in drosophila eye: I dynamics. PLoS One, 4, e4307.
Contents lists available at ScienceDirect
Neural Networks journal homepage: www.elsevier.com/locate/neunet
The transfer function of neuron spike Igor Palmieri a,∗ , Luiz H.A. Monteiro b,a , Maria D. Miranda a a
Departamento de Engenharia de Telecomunicações e Controle, Escola Politécnica da Universidade de São Paulo, Av. Prof. Luciano Gualberto, travessa 3, n. 158, 05508-900, São Paulo, SP, Brazil b Escola de Engenharia da Universidade Presbiteriana Mackenzie, Pós-graduação em Engenharia Elétrica, Rua da Consolação, n. 896, 01302-907, São Paulo, SP, Brazil
article
info
Article history: Received 20 December 2014 Received in revised form 9 March 2015 Accepted 8 April 2015 Available online 22 April 2015 Keywords: Action potential Neuronal spike Parametric modeling Transfer function
abstract The mathematical modeling of neuronal signals is a relevant problem in neuroscience. The complexity of the neuron behavior, however, makes this problem a particularly difficult task. Here, we propose a discrete-time linear time-invariant (LTI) model with a rational function in order to represent the neuronal spike detected by an electrode located in the surroundings of the nerve cell. The model is presented as a cascade association of two subsystems: one that generates an action potential from an input stimulus, and one that represents the medium between the cell and the electrode. The suggested approach employs system identification and signal processing concepts, and is dissociated from any considerations about the biophysical processes of the neuronal cell, providing a low-complexity alternative to model the neuronal spike. The model is validated by using in vivo experimental readings of intracellular and extracellular signals. A computational simulation of the model is presented in order to assess its proximity to the neuronal signal and to observe the variability of the estimated parameters. The implications of the results are discussed in the context of spike sorting. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction The mathematical modeling of neuronal signals is one of the most relevant problems in neuroscience. This question has been driven by technological advances that allowed the experimental acquisition of neuronal signal in several situations, establishing the basis for many quantitative models found in literature (Koch & Segev, 1998; Sterratt, Graham, Gillies, & Willshaw, 2011). Such models have been widely employed to solve computational tasks in a computer (Gao, Yang, Cai, & Liu, 2012; Menezes & Monteiro, 2011) and to understand the functioning of nervous systems (Monteiro, Bussab, & Chaui-Berlinck, 2002; Steuber & Jaeger, 2013). One of the first models was proposed by Huxley and Hodgkin in the 1950s, using a system of continuous-time equations based on experiments with giant squid cells (Hodgkin & Huxley, 1952). This model served as inspiration for subsequent proposals, mainly focused on the reproduction of the action potential dynamics, as in FitzHugh (1961). Discrete-time versions (Ibarz, Casado, & Sanjuán, 2011), as the McCulloch–Pitts model, can be considered simplifications to facilitate simulations in digital environments.
∗
Corresponding author. E-mail addresses: [email protected] (I. Palmieri), [email protected] (L.H.A. Monteiro), [email protected] (M.D. Miranda). http://dx.doi.org/10.1016/j.neunet.2015.04.003 0893-6080/© 2015 Elsevier Ltd. All rights reserved.
There are many biophysical and chemical factors involved in the generation of neuronal signals; hence, to model such signals is a particularly difficult task. Moreover, the neuron behavior is highly influenced by the connections with surrounding nerve cells, forming a network which itself poses a great challenge for mathematical representation (Dayan & Abbott, 2005). Most of the referred models are designed with the main intent of representing the neuronal dynamics and the interaction with external agents. However, in several situations, a parametrization associated solely with the signal waveform is of great interest. Extra and intracellular action potential recordings have been used to tune the parameter values of continuous-time nonlinear models, which explicitly consider the physical properties of the neuronal membrane. This tuning is commonly performed by using parameter estimation techniques (Gold, Henze, Koch, & Buzsáki, 2006; Tabak, Murphey, & Moore, 2000; Toth, Kostuk, Meliza, Margoliash, & Abarbanel, 2011; Wang et al., 2014; Willms, Baro, Harris-Warrick, & Guckenheimer, 1999). However, the number of parameters to be matched can be as great as 100 (Gold, Henze, & Koch, 2007). In fact, the complexity of these models can impair computational simulations and analytical studies involving many neurons. There are also discrete-time nonlinear models that can generate a signal which is qualitatively similar to neuronal recordings (Aihara, Takabe, & Toyoda, 1990; GirardiSchappo, Tragtenberg, & Kinouchi, 2013; Hsu, Gobovic, Zaghloul,
90
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
& Szu, 1996; Mesbah, Moghtadaei, Golpayegani, & Towhidkhah, 2014). However, these models neither take into account details of the neurophysiological processes nor experimental data. Their parameter values are empirically set. The parametrization of intracranial signals is of great interest in the context of the identification problem. The study of efficient techniques to isolate the firings of each active cell from an extracellular mixture containing sequences of spikes from several neurons is a problem extensively discussed in Neuroscience, formally known as spike sorting. Usually these techniques are implemented in three stages (Gibson, Judy, & Markovic, 2012; Lewicki, 1998; Quiroga, Nádasdy, & Ben-Shaul, 2004; Rutishauser, Schuman, & Mamelak, 2006): from raw electrical potentials recorded by using arrays of electrodes, spikes are detected, then parameterized, and finally sorted, attributing every single spike observed to a particular neuron. In this context, a relevant question is how the parameter sets describe effectively the neuronal activities over time. A simplified model of how the signal that arrives at the electrodes is formed from a stimulus and its evolution over time can possibly help to elucidate this question. The main goal of this paper is to propose a discrete-time linear time-invariant (LTI) model with a rational function, to represent the neuronal spike observed by an electrode located in the surroundings of the nerve cell. This approach provides a low-complexity alternative to the neuronal spike model, employing system identification and signal processing concepts. It is dissociated from any considerations about the biophysical processes of the neuronal cell, and it is not directly derived from any traditional neuronal models. The model is validated by using in vivo experimental readings of intracellular and extracellular signals. Although the proposed model does not fully comprise the cellular dynamics, it allows a straightforward simulation of the neuronal activity by a pair of linear equations, once its parameters are obtained. The remainder of this paper is organized into five sections. In Section 2, the experimental methods involved in the data acquisition of neuronal electrical signals are briefly described. The proposed model is introduced in Section 3, in the form of a cascade association of two subsystems: one that generates an action potential from an input stimulus, and one that represents the medium between the cell and the electrode. The parameters for these subsystems are estimated from experimental data in Section 4. In Section 5, a computational simulation of the model is presented in order to assess its proximity to the neuronal signal and the variability of the estimated parameters. The implications of the results are discussed in Section 6, considering the limitations of the proposed approach. 2. Experimental methods The estimation of the proposed model parameters is performed by using in vivo experimental data recordings of neuronal signals generated by mice, provided by CRCNS (Collaborative Research in Computational Neuroscience) (Henze, Harris et al., 2000). For each experiment in the database, there are time series of action potentials obtained directly from the cell body through glass micropipettes, and of the electric potential observed simultaneously in the external environment, through electrodes placed in an area near the chosen cell. These external electrodes are mounted in arrays of four to six uniformly spaced elements, providing spatial diversity in the recordings of the electric potential outside the cell. All the recordings were obtained from the region of the hippocampus; a great part of the pulses corresponds to the pyramidal cell activity observed in a layer known as CA1. The time series were digitized at 12 or 16 bits, and sampled at frequency
Fig. 1. (a) Segment of the original acquired extracellular signal; (b) the same segment after filtering to remove low frequency components.
of 10, 20, 25 or 50 kHz, depending on the experiment. The values were then recorded in binary files, and made publicly available by CRCNS at their website (Henze, Harris et al., 2000). According to the coordinators of the experiment, the extracellular signal was first amplified 1000–8000 times, and filtered at the source by a low pass filter with cutoff frequency of 3 kHz (Henze, Borhegyi et al., 2000). Before using the digitalized data in the modeling process, however, we also employed a filtering process in each sequence. In the extracellular sequence, we applied an elliptical-type high-pass filter with cutoff frequency of 300 Hz. The purpose of this filtering is to remove low frequency components, mainly associated with the LFP (Local Field Potential). Although it also contains some information about the neuronal activity (Scherberger, Jarvis, & Andersen, 2005), this fact is not addressed in the proposed parametrization. Fig. 1 shows a fragment of the extracellular signal before and after this stage. The intracellular sequence, composed almost exclusively by the action potential of the neuron, is also passed through a high-pass filter, but with cutoff frequency of 50 Hz. This filtering only aims to remove low frequency components, especially in zero frequency, to ensure the level of the resting potential near zero. After the filtering process, we selected a fixed interval of 25 s within one of the available recordings, containing 86 action potentials of the same neuron. We noted that the corresponding intracellular series contains solely these action potentials, as it is obtained directly in the membrane of a single neuron, while the extracellular series contains the spikes of all neurons in the region during the same time interval. Centered around the exact sample of the maximum value of each signal, we defined a window of 14 ms (N = 140 samples for a sampling rate of 10 kHz). In the remainder of this paper, the intracellular and extracellular signals inside each time window are denoted as {s(ℓ)} and {y(ℓ)}, respectively, and the estimate of these signals are denoted as {ˆs(ℓ)} and {ˆy(ℓ)} with ℓ = 0, . . . , N − 1. 3. The proposed model The main goal of parametric modeling is to build a mathematical model able to represent relevant characteristics of a given signal by a reduced set of parameters. Here we consider the specific case of modeling the discrete-time extracellular signal y(n) whose sample set is obtained through direct observation inside a fixed time window. Fig. 2 shows the simplified block diagram of the proposed model, in which x(n) is the input and the output yˆ (n) is the estimate of y(n). We assume that this signal can be represented by a LTI
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
91
Fig. 2. Simplified block diagram of the proposed model: F (z ) is the all-pole model, G(z ) is the moving average model, yˆ (n) is the output signal, and η(n) is an additive noise include to yield a signal similar to the recordings.
model with transfer function H (z ). This function is represented with two subsystems F (z ) and G(z ) associated in series, so that q
gm
b k z −k
k=0 p
H (z ) , F (z ) G(z ) =
1−
(1) ak z − k
k=1
in which gm is a gain factor, and p and q represent the orders of the polynomials at the denominator and at the numerator, respectively. The parametrization process corresponds to estimate H (z ), or the values of the coefficients ak and bk (and also the gain gm ), that together define a model capable of yielding an approximation yˆ (n) of the original signal y(n), given an input x(n). Essentially, when a punctual stimulus is above a threshold, the neuron fires. If the stimulus received by the neuron is below such a threshold, it can be considered null; if it is above, its value is taken as equal to one. This assumption is common in artificial neural networks, where the neuron output is usually a binary variable; typically, 0 or 1 (Haykin, 2009). As a premise, in our approach, it is considered that y(n) is basically a deterministic signal, and the input of the model is a Kronecker delta, or discrete unit impulse (Hayes, 1996), time shifted to the exact moment when the neuron fired. This unit impulse x(n) = δ(n − ∆) =
1, 0,
n=0 n ̸= 0
(2)
is essentially a punctual stimulus on the nerve cell enough to make it fire at n = ∆. So, in our approach, the estimate of the observed extracellular spike output yˆ (n) has the same role as the impulse response of the system. Hence, H (z ) can be considered the transfer function of this system. The subsystem F (z ) represents the model associated with the action potential produced in the cell membrane, and thus it represents characteristics inherent to each neuron. We assume that the input–output relation is an autoregressive all-pole model, containing the denominator of the transfer function, scaled by a gain factor gm , that is, gm
F (z ) = 1−
p
.
(3)
ak z − k
k=1
Then the output sˆ(n) generated by this model is calculated recursively by the relation sˆ(n) =
p k=1
ak sˆ(n − k) + gm x(n) = sˆ Tp (n − 1) ap + gm x(n)
(4)
in which ap = [a1 a2 . . . ap ]T and sˆ p (n − 1) = [ˆs(n − 1) sˆ(n − 2) . . . sˆ(n − p)]T are the coefficients and the input column vectors respectively, and T denotes the vector transpose operator. In broad terms, the output sˆ(n) should be a good estimate of the original intracellular signal s(n). If the signal s(n) could be entirely represented by an all-pole system, so a proper choice of parameters would result in sˆ(n) = s(n). The subsystem G(z ) represents the medium between the cell and the electrode, and it acts as a communication channel through which the generated signal propagates. We assume that the input–output relation is a moving average model, that is, G(z ) =
q
b k z −k .
(5)
k=0
Then, its output is calculated by the relation yˆ (n) =
q
bk sˆ(n − k) = sˆ Tq+1 (n) bq+1
(6)
k=0
in which bq+1 = [b0 b1 . . . bq ]T and sˆ q+1 (n) = [ˆs(n) sˆ(n − 1) . . . sˆ(n − q)]T are the coefficients and the input column vectors, respectively. It is desired that the output of the complete system, formed by the concatenation of F (z ) with G(z ), provides an estimate of the signal y(n), recorded from the outside medium of the cell, with characteristics that depend on the neuron itself and the environment shared by the cells in that specific brain region. The output of the system formed by the cascade of the autoregressive and moving average models comprises the deterministic part of the signal that we are trying to model. However, observations of actual extracellular signal indicate the presence of some level of noise η(n) due to the measurement process. This effect may be represented by adding a random signal at the output of the function of the model. This noise can also take in account the variability in the underlying neural circuits. 4. Parameter estimation of the transfer function There are several approaches that can be used to calculate the coefficients of the model in order to approximate the actual signal in each block of the unknown system (Ljung, 1999). This work considers Prony’s method (Hayes, 1996) to produce the estimates of the autoregressive coefficients, associated with the model F (z ) of intracellular signal. In the case of the model G(z ), the optimal coefficients are estimated by using a system identification configuration whose input is the intracellular signal s(n), and the output is the extracellular signal y(n). In both cases we apply the Least Squares Error (LSE) criterion. For each modeled action
92
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
Fig. 3. Average intracellular signal and limits defined by the standard deviation.
Fig. 4. Average extracellular signal and limits defined by the standard deviation.
potential and spike, only a finite number N of samples are available. In our approach for estimated all the coefficients of F (z ) and G(z ), we consider the autocorrelation method which makes the assumptions that the input data prior to n = 0 and the data after n = N − 1 are zero (Oppenheim & Schafer, 2009).
this in mind, we define the estimation error between the desired signal y(n) and model output yˆ (n) as ey (n) = y(n) − yˆ (n)
(10)
in which yˆ (n) = sˆ Tq+1 (n) bq+1 . The solution for the best transversal filter, according to LSE criterion, is given by the equation
4.1. The parameters of F (z )
b q +1 =
N +q
−1 sq+1 (ℓ)
sTq+1
(ℓ)
N +q
sq+1 (ℓ)y(ℓ).
(11)
Prony’s method is a way of estimating the coefficients of an autoregressive and moving-average model, especially useful when applied to the case of the deterministic all-pole model (Hayes, 1996). The coefficients are obtained by solving the problem of forward prediction of the desired signal. The forward prediction error is defined as
Here, the inverse of the autocorrelation matrix is estimated by using the experimental readings of the intracellular signal {s(ℓ)}, as in Eq. (9), but the cross-correlation vector is estimated by using the intracellular {s(ℓ)} and extracellular {y(ℓ)} signals.
es (n) = s(n) − sˆ(n)
5. Results and simulation
(7)
in which sˆ(n) = (n − 1)ap . The optimal coefficients are obtained in order to minimize the energy of es (n). A compact way to represent this criterion is sTp
ap = arg min ap
N −1
|es (n)|
2
(8)
n =0
where ap is the vector of optimal coefficients. The solution that satisfies (8), according to the LSE criterion, must also satisfy the relation
ap =
−1
N +p−1
sp (ℓ) (ℓ) sTp
N +p−1
(9)
ℓ=0
ℓ=0
sp (ℓ)s(ℓ + 1) .
(I )
(II )
It is worth to note that the inverse of the autocorrelation matrix (I ) and the crosscorrelation vector (II ) are estimated by using in vivo experimental readings of intracellular signals {s(ℓ)}. The value of gm is calculated in a later step, as a gain factor to equalize the energy of the observed signal and the estimated signal. The last task of intracellular signal modeling is to define exactly the input x(n). For this end, an integer delay value ∆m is determined empirically, to be applied on the input signal, such that x(n) = δ(n − ∆m ) acts as a trigger in the same instant that the real neuron fired. 4.2. The parameters of G(z ) Given the input s(n) and the output y(n) of an unknown system, the system identification configuration allows us to obtain the coefficients bk of the moving average model (Haykin, 1996). With
ℓ=0
ℓ=0
In the following, we present the selected firings contained in the filtered signals from the database. Then, we show how the orders of the subsystems F (z ) and G(z ) are estimated. Finally, all the parameters values are calculated. With these coefficients, we obtain the estimates sˆ(n) and yˆ (n), which are compared with the original signals s(n) and y(n). 5.1. In vivo experimental readings of intracellular and extracellular signals We selected a 25 s segment of the intra and extracellular signals filtered from the original time series of the database, containing a sequence of 86 action potentials. Each action potential is represented by a window with N = 140 samples (14 ms). Despite the fact that the firing sequence has been generated by the same neuron, the waveform of the action potential and the spike do not have a fixed shape over time. Figs. 3 and 4 illustrate the mean values of intracellular and extracellular signals inside the 14 ms window, considering all the firing instants in the chosen segment. Observe that there is a noticeable variability in the waveform of both signals, and also an apparent level of noise in the extracellular sequence. 5.2. The estimate of the model order For the subsystem F (z ) and G(z ), it is necessary to define a number of coefficients that allow the model to reproduce the signal of interest with reasonable assertiveness. In other words, we must find the integers p and q that suffice to build a suitable model. This task is directly linked with the complexity of the optimal coefficients estimations in Eqs. (9) and (11).
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
93
Fig. 5. Partial autocorrelation function of the intracellular signal, with a confidence interval based on the standard error (SE) to help the definition of the model order. Fig. 7. Model output sˆ(n) vs. acquired intracellular signal s(n).
Fig. 6. Model order vs. normalized (by the energy of the real signal) mean squared error. Fig. 8. Model output yˆ (n) vs. acquired extracellular signal y(n).
In the all-pole model case, a good indication of the number of coefficients of the F (z ) is given by the partial autocorrelation function of the signal being modeled, while its value is far from zero (Ljung, 1999). This measure, however, does not make sense to the moving average extracellular model G(z ). As an alternative, in both subsystems, it is also possible to calculate the normalized mean squared error between the model output and the observed signal as function of the model order (Oppenheim & Schafer, 2009). When this error approaches a constant level, there is an indication that the addition of coefficients is no longer relevant. Fig. 5 shows the average partial autocorrelation function calculated with the intracellular time series, and the 95% confidence interval of the standard error (Box, Jenkins, & Reinsel, 2008). It is also possible to simulate both models with a different number of coefficients and observe their residual error. Fig. 6 exhibits the normalized mean squared error of intracellular and extracellular signals (given by Eqs. (7) and (10), respectively) of a single neuron, compared to the model order. These curves together suggest that we can use p = 6 and q = 3 coefficients for the intra and extracellular models, respectively. 5.3. The estimate of the intra and extracellular signals The samples of each window are then used to estimate the autocorrelation and the crosscorrelation functions, which will solve Eqs. (9) and (11). Each window therefore corresponds to a single set of model parameters, for both the intracellular and the extracellular models. We can better observe the result of the modeling process by analyzing the time interval of a single isolated action potential and its corresponding spike. Fig. 7 illustrates the output produced by the all-pole intracellular model along with the signal used to estimate the coefficients. We can see that there are some differences in the generated signal, especially in the repolarization phase.
The moving average extracellular subsystem is also simulated within the same time frame. In Fig. 8, the output of the complete model is shown, compared to the observed extracellular spike. Note that there is an influence of the experimental noise in this signal, which accounts for the main differences from the model output. 5.4. Distributions of the coefficients The variability of the waveform in the intra and extracellular signals is reflected in the changes of the estimates of the F (z ) and G(z ) coefficients, respectively. Figs. 9 and 10 present the normalized histogram (such that the sum of the bars is equal to one) of the coefficients calculated for each subsystem using the chosen sequence of firings. It is possible to observe the behavior of the model coefficients for multiple firings in the same sequence by looking at the statistical distribution of their values. Another important result is that the average coefficients of the extracellular model, seen in Fig. 10, indicate a transfer function of a differentiator, as expected given the capacitive characteristics of the external medium through which the signal propagates (Gibson et al., 2012). By inserting these average coefficients values in Eq. (4), with the average gain, we obtain the following relation to calculate the output for the observed neuron: sˆ(n) = sˆ T6 (n − 1) a6 + 14.622 δ(n − ∆)
(12)
in which a6 = [2.436 − 2.430 1.456 − 0.714 0.267 − 0.045]T . These coefficients form a system whose poles are inside the unitary circle, and consequently this is a bounded-input, boundedoutput (BIBO) stable system. Considering that the neuron fired at the instant n = ∆, the resulting values for the output will start at zero and then follow the sequence sˆ(∆ − 1) = 0; sˆ(∆) = 14.622; sˆ(∆ + 1) = 35.619; sˆ(∆ + 2) = 51.237, and so on, decreasing asymptotically to the initial signal level.
94
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
types of neurons (Bean, 2007). Thus, although Eqs. (12) and (13) were obtained from a particular database, the method here proposed can be used for modeling spikes with distinct waveforms. This is relevant because neural communication is based not only on frequency and phase of action potentials, but also in their amplitude (Bean, 2007; de Polavieja et al., 2005; Zheng et al., 2009). 6. Discussion
Fig. 9. Distribution of the estimated coefficients of the intracellular model, with average values indicated.
The model proposed in this paper has a high level of abstraction regarding any physiological process and presents a considerable level of simplification, being comparable to most of the discretetime models found in the literature (Dayan & Abbott, 2005). Its main feature is to employ a linear equation for representing the neuronal signal during the firing of the cell, instead of a set of non-linear equations. The output signal of the model is readily calculated once its parameters are estimated. Usually, neuronal models try to reproduce or mimic the dynamics involved in the generation of an action potential, including features like fire rate or the cell’s membrane conductances (Ibarz et al., 2011; Izhikevich, 2007). In the proposed model, the goal is to describe mainly its waveform. The generation of the action potential and the resulting spike are triggered by the input signal, containing a unit impulse time shifted to the exact moment of the firing. Despite showing minor deviations from the original signal, the main significance of this work is the structure for the neuronal signal model: autoregressive in the intracellular and moving average in the extracellular. It should also be noted that the parameters of the transfer function H (z ) vary between each firing, reflecting the changes in the signal waveform. These variations are observed in both intra and extracellular signals, even though they present different signal-to-noise ratio. Expected developments
Fig. 10. Distribution of the estimated coefficients of the extracellular model, with average values indicated.
An analogous relation can be calculated for the extracellular model, using the average coefficients in the Eq. (6): yˆ (n) = sˆ T4 (n) b4
(13)
in which b4 = [−0.454 0.452 − 0.275 0.228] . This signal also departs from zero, and the input remains zero till the neuron fires. The output sequence will follow the values yˆ (∆ − 1) = 0; yˆ (∆) = −6.638; yˆ (∆ + 1) = −9.562; yˆ (∆ + 2) = −11.183, also decreasing back to the same initial values. The dispersion around their mean values is a consequence of the signal variability, as displayed in Figs. 3 and 4. Such changes in the shape of action potentials can be due to differences in ionic conductances (Djouhri, Bleazard, & Lawson, 1998); changes on the levels of neurohormones, neuromodulators, neurotransmitters (Rodríguez-Molina et al., 2014); diversity on the stimulus history (de Polavieja, Harsch, Kleppe, Robinson, & Juusola, 2005); distinct T
In practical applications, an electrode employed in neural signal acquisition registers the combined activity of a set of neighbor neurons. The goal of spike sorting processes is to distinguish the neurons generating the activity detected by every electrode. This task includes stages of clustering and classification of the detected spikes, and is highly based on a parametric representation of the extracellular signal (Gibson et al., 2012; Lewicki, 1998; Quiroga et al., 2004; Rutishauser et al., 2006; Sarinho Filho, Eisencraft, Suyama, Fonoff, & Miranda, 2011). The results obtained with the proposed approach suggest further investigation of the connection between the polynomial H (z ) and the parametrization methods of spike sorting. However, simulations showed that the changes in the spike waveform are directly reflected in the model coefficients, even when the presence of noise is ignored. Consequently, it should be expected that any method of spike identification based on the extracellular signals from experimental data will also subject to the effects of variations in waveform parameters. Due to changes in the waveform of the spike, errors in parametrization and classification phases of any spike sorting process are hard to overcome. This suggests that the dynamics of the neuronal activity can be better described if the decision about the neuron activation is observed over time consecutively with the detection stage and not in a posterior stage like is commonly done in spike sorting. In such context, the use of adaptive techniques combined with the detection and the parameter estimation stages is of high interest, as addressed in Miranda and Gerken (1997) and Miranda, Silva, and Nascimento (2008). A deeper investigation should explore the role of the model parameters in the process of neuron identification.
I. Palmieri et al. / Neural Networks 68 (2015) 89–95
The model can also be further expanded by considering simultaneous readings of different extracellular electrodes, resulting in a system containing multiple output channels. Since usually the electrodes are located in different points of the external medium, the resulting model should reflect the spatial diversity of the extracellular signal. A more general structure should also include more than one input, each one associated with the firings of an active cell in a region containing several neurons. However, this model would require observations of the intracellular signal of other cells. Due to the availability of the experimental data, consisting of only one intracellular time series, the simulations in this work were so far restricted to the single input case. Acknowledgment L.H.A. Monteiro is partially supported by National Council for Scientific and Technological Development (CNPq) under the grant #305827/2014-6. References Aihara, K., Takabe, T., & Toyoda, M. (1990). Chaotic neural networks. Physics Letters A, 144, 333–340. Bean, B. P. (2007). The action potential in mammalian central neurons. Nature Reviews Neuroscience, 8, 451–465. Box, G. E. P., Jenkins, G. M., & Reinsel, G. C. (2008). Time series analysis: forecasting and control (4th ed.). Hoboken, NJ: John Wiley & Sons. Dayan, P., & Abbott, L. F. (2005). Theoretical neuroscience: computational and mathematical modeling of neural systems. Cambridge, MA: MIT Press. de Polavieja, G. G., Harsch, A., Kleppe, I., Robinson, H. P. C., & Juusola, M. (2005). Stimulus history reliably shapes action potential waveforms of cortical neurons. The Journal of Neuroscience, 25, 5657–5665. Djouhri, L., Bleazard, L., & Lawson, S. N. (1998). Association of somatic action potential shape with sensory receptive properties in guinea-pig dorsal root ganglion neurones. The Journal of Physiology, 513, 857–872. FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1, 445–466. Gao, D., Yang, Z., Cai, C., & Liu, F. (2012). Performance evaluation of multilayer perceptrons for discriminating and quantifying multiple kinds of odors with an electronic nose. Neural Networks, 33, 204–215. Gibson, S., Judy, J. W., & Markovic, D. (2012). Spike sorting: the first step in decoding the brain. IEEE Signal Processing Magazine, 29, 124–143. Girardi-Schappo, M., Tragtenberg, M. H. R., & Kinouchi, O. (2013). A brief history of excitable map-based neurons and neural networks. Journal of Neuroscience Methods, 220, 116–130. Gold, C., Henze, D. A., & Koch, C. (2007). Using extracellular action potential recordings to constrain compartmental models. Journal of Computational Neuroscience, 23, 39–58. Gold, C., Henze, D. A., Koch, C., & Buzsáki, G. (2006). On the origin of the extracellular action potential waveform: a modeling study. Journal of Neurophysiology, 95, 3113–3128. Hayes, M. H. (1996). Statistical digital signal processing and modeling. New York, NY: John Wiley & Sons. Haykin, S. (1996). Adaptive filter theory (3rd ed.). Upper Saddle River, NJ: Prentice Hall. Haykin, S. (2009). Neural networks: a comprehensive foundation. Upper Saddle River, NJ: Prentice Hall. Henze, D. A., Borhegyi, Z., Csicsvari, J., Mamiya, A., Harris, K. D., & Buzsáki, G. (2000). Intracellular features predicted by extracellular recordings in the hippocampus in vivo. Journal of Neurophysiology, 84, 390–400. Henze, D.A., Harris, K.D., Borhegyi, Z., Csicsvari, J., Mamiya, A., & Hirase, H. et al. 2000. Simultaneous intracellular and extracellular recordings from hippocampus region CA1 of anesthetized rats. CRCNS.org. URL: http://dx.doi.org/10.6080/K02Z13FP.
95
Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117, 500–544. Hsu, C. C., Gobovic, D., Zaghloul, M. E., & Szu, H. H. (1996). Chaotic neuron models and their VLSI circuit implementations. IEEE Transactions on Neural Networks, 7, 1339–1350. Ibarz, B., Casado, J. M., & Sanjuán, M. A. F. (2011). Map-based models in neuronal dynamics. Physics Reports, 501, 1–74. Izhikevich, E. M. (2007). Dynamical systems in neuroscience. Cambridge, MA: MIT Press. Koch, C., & Segev, I. (Eds.) (1998). Methods in neuronal modeling: from ions to networks (2nd ed.). Cambridge, MA: MIT Press. Lewicki, M. S. (1998). A review of methods for spike sorting: the detection and classification of neural action potentials. Network: Computation in Neural Systems, 9, 53–78. Ljung, L. (1999). System identification: theory for the user (2nd ed.). Upper Saddle River, NJ: Prentice Hall. Menezes, R. A., & Monteiro, L. H. A. (2011). Synaptic compensation on Hopfield network: implications for memory rehabilitation. Neural Computing and Applications, 20, 753–757. Mesbah, S., Moghtadaei, M., Golpayegani, M. R. H., & Towhidkhah, F. (2014). Onedimensional map-based neuron model: a logistic modification. Chaos, Solitons and Fractals, 65, 20–29. Miranda, M. D., & Gerken, M. (1997). A hybrid least squares QR-lattice algorithm using a priori errors. IEEE Transactions on Signal Processing, 45, 2900–2911. Miranda, M. D., Silva, M. T. M., & Nascimento, V. H. (2008). Avoiding divergence in the Shalvi–Weinstein algorithm. IEEE Transactions on Signal Processing, 56, 5403–5413. Monteiro, L. H. A., Bussab, M. A., & Chaui-Berlinck, J. G. (2002). Analytical results on a Wilson–Cowan neuronal network modified model. Journal of Theoretical Biology, 219, 83–91. Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-time signal processing (3rd ed.). Upper Saddle River, NJ: Prentice Hall. Quiroga, R. Q., Nádasdy, Z., & Ben-Shaul, Y. (2004). Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering. Neural Computation, 16, 1661–1687. Rodríguez-Molina, V., Patiño, J., Vargas, Y., Sánchez-Jaramillo, E., Joseph-Bravo, P., & Charli, J. L. (2014). TRH regulates action potential shape in cerebral cortex pyramidal neurons. Brain Research, 1571, 1–11. Rutishauser, U., Schuman, E. M., & Mamelak, A. N. (2006). Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo. Journal of Neuroscience Methods, 154, 204–224. Sarinho Filho, J. N. S., Eisencraft, M., Suyama, R., Fonoff, E. T., & Miranda, M. D. (2011). The use of least squares lattice algorithm in the parameterization and sorting of action potentials signals. In Bioelectronics, biomedical, and bioinspired systems V and nanotechnology V . 80680Q. Scherberger, H., Jarvis, M. R., & Andersen, R. A. (2005). Cortical local field potential encodes movement intentions in the posterior parietal cortex. Neuron, 46, 347–354. Sterratt, D., Graham, B., Gillies, A., & Willshaw, D. (2011). Principles of computational modelling in neuroscience. Cambridge, UK: Cambridge University Press. Steuber, V., & Jaeger, D. (2013). Modeling the generation of output by the cerebellar nuclei. Neural Networks, 47, 112–119. Tabak, J., Murphey, C. R., & Moore, L. (2000). Parameter estimation methods for single neuron models. Journal of Computational Neuroscience, 9, 215–236. Toth, B. A., Kostuk, M., Meliza, C. D., Margoliash, D., & Abarbanel, H. D. I. (2011). Dynamical estimation of neuron and network properties I: variational methods. Biological Cybernetics, 105, 217–237. Wang, R., Wang, J., Deng, B., Liu, C., Wei, X., Tsang, K. M., et al. (2014). A combined method to estimate parameters of the thalamocortical model from a heavily noise-corrupted time series of action potential. Chaos, 24, 013128. Willms, A. R., Baro, D. J., Harris-Warrick, R. M., & Guckenheimer, J. (1999). An improved parameter estimation method for Hodgkin–Huxley models. Journal of Computational Neuroscience, 6, 145–168. Zheng, L., Nikolaev, A., Wardill, T. J., O’Kane, C. J., de Polavieja, G. G., & Juusola, M. (2009). Network adaptation improves temporal representation of naturalistic stimuli in drosophila eye: I dynamics. PLoS One, 4, e4307.
