An Algorithm for Tracking the Position and Velocity of Multiple Neuronal Signals Using Implantable Microelectrodes In Vivo

Abstract

Increasingly complex multi-electrode arrays for the study of neurons both in vitro and in vivo have been developed with the aim of tracking the conduction of neural action potentials across a complex interconnected network. This is usually performed through the use of electrodes to record from single or small groups of microelectrodes, and using only one electrode to monitor an action potential at any given time. More complex high-density electrode structures (with thousands of electrodes or more) capable of tracking action potential propagation have been developed but are not widely available. We have developed an algorithm taking data from clusters of electrodes positioned such that action potentials are detected by multiple sites, and using this to detect the location and velocity of action potentials from multiple neurons. The system has been tested by analyzing recordings from probes implanted into the locust nervous system, where recorded positions and velocities correlate well with the known physical form of the nerve.

1. Introduction

Since the first development of action potential measurement in the 1950s [1], instruments have been developed to enable the monitoring of single action potentials (APs). Whilst original recording devices comprised thin wires, since the 1970s, the predominant approach to neural recording electrodes has been through the use of micromachined silicon electrodes constructed using techniques developed in the semiconductor industry. As manufacturing techniques have improved, it has been possible to increase the number of electrodes on each device in order to enable recording from multiple neurons simultaneously, either in in vivo or in vitro settings [2,3,4,5]. When recording from complex structures containing many neurons, it is important to be able to discriminate between different neural sources. Where single electrodes are used in multi-neuron settings, this is commonly performed by classifying the sources according to amplitude to determine proximity, though this is only effective for AP sources adjacent to the electrode; as distance increases, so does the likelihood of multiple sources being at similar distance from the electrode in any given direction; such single electrode systems (or multi-electrode systems with inter-electrode spacings of hundreds of microns) are unable to give spatial information about these sources, such as location or velocity. Pickard et al. demonstrated that a single implanted multi-electrode probe can detect many action potentials by comparing the relative magnitudes of the potential “spike” [6]. This approach has been used by other groups to discriminate between APs, but the typical wide spacing between microelectrodes of implantable devices (typically multiple hundreds of microns) means that AP spikes can be detected by only one electrode at a given point in time (e.g., [7]).

Where spatial information such as AP location, detection and velocity are required, multiple electrodes must be used. A common approach to this is to track APs using large multiple-electrode arrays (MEAs), such that individual APs pass near many electrodes and can consequently be identified, tracked and measured. However, such MEAs typically use high levels of redundancy, with many electrodes (sometimes extending to hundreds) required to track each AP, and each electrode recording from a small number of sources [8]. Whilst in vitro MEAs can be manufactured to fit the size of a brain slice or cultured neural structure, in vivo devices are more limited in order to maximize proximity to neuron sources without causing significant damage to the neural tissue. Consequently, for in vivo electrodes, a more efficient approach to detection of APs is required, using smaller numbers of electrodes in close proximity to one another and use triangulation techniques in order to identify multiple AP sources. This has been used to measure the position of multiple sources in two dimensions using parallel metal wires [9], though the nature of the very long electrodes placed parallel to the direction of neural processes means the approach is only able to classify forces radially, and relies on the electrode direction being completely parallel to the neural sources. Another approach developed by Hughes et al. [10] used point electrodes to measure action potential velocity by measuring the time taken to travel between electrodes placed along the nerve.

In this paper, we present an approach to detecting action potentials from multi-unit recordings that allow us to track the axon source across a nerve with as few as three electrodes by triangulating the signals from multiple electrodes, and comparing their timings and amplitudes. In order to evaluate the usefulness of the algorithm, we have validated the model by applying it to signals recorded using a microelectrode device implanted into a locust ventral nerve cord, a complex network of non-myelinated nervous systems with defined physical dimensions and with action potential travelling along a defined afferent-efferent direction of travel. The result demonstrated a clear correlation between the analyzed results and the nerve geometry. This has the potential to provide a new approach to multi-unit neural recordings, simplifying the production of electrode arrays for the analysis of very large numbers of neurons without the requirement for commensurately large numbers of electrodes.

2. Detection Algorithm

We aimed to develop a method of AP tracking that minimizes the number of electrodes required for neural monitoring, by using three electrodes arranged on a penetrating probe with a maximum distance between the furthest electrodes of less than 100 µm. These were used not only to identify individual APs, but also to calculate their position, direction and velocity. The analysis approach is derived from a model of the movement of an AP past multiple electrodes.

In order to explain the analysis method, we will describe the model on which it is based, which can be seen in Figure 1. At its simplest, we can model an AP as a charge Q travelling in a linear axon parallel to the probe surface. From Coulomb’s law, the potential caused by a charge at a point a will be inversely proportional to the distance between that point and the charge. At the scale of the probe, a nerve can reasonably be considered linear and the position of the probe can be made as to satisfy the condition of parallelism.

As our electrodes are in planar arrays, we can define them purely in the x–y plane. We can simplify the electrodes by considering them as a recording point with position ${\overrightarrow{M}}_i=x_i.{\overrightarrow{u}}_x+y_i.{\overrightarrow{u}}_y$. The potential at each electrode at time t depends on the distance D_i between the charge and the electrode, which can be expressed as a function of the position s of the charge along the axon and on the minimum distance between the electrode and the axon. Taking into account the AP velocity $\overrightarrow{v}=v_x.{\overrightarrow{u}}_x+v_y.{\overrightarrow{u}}_y$, one obtains the following:

$${\overrightarrow{\mathit{D}}}_{\mathit{i}}=-{\overrightarrow{\mathit{M}}}_{\mathit{i}}+{\overrightarrow{\mathit{s}}}_{\mathit{o}}+\mathit{t}\cdot \overrightarrow{\mathit{v}}$$

(1)

where

${\overrightarrow{s}}_0=x_0.{\overrightarrow{u}}_x+y_0.{\overrightarrow{u}}_y+z_0.{\overrightarrow{u}}_z$ is a reference point on the axon defined by the position of the AP at time 0.

The potential at the electrode reaches its maximum value when the charge is at the point of the axon closest to the probe. We call this point ${\overrightarrow{s}}_{m,i}$ with coordinates $x_{m,i}, y_{m,i}, z_{m,i}$, where m corresponds to ‘minimum’ and i describes reference number of the electrode. Since the probe is parallel to the axon, one can write:

$${\overrightarrow{\mathit{M}}}_{\mathit{i}}·{\overrightarrow{\mathit{s}}}_{\mathit{m},\mathit{i}}=\mathbf{0}=\left({\mathit{x}}_{\mathit{m},\mathit{i}}-{\mathit{x}}_{\mathit{i}}\right)·\mathit{t}·{\mathit{v}}_{\mathit{x}}+\left({\mathit{y}}_{\mathit{m},\mathit{i}}-{\mathit{y}}_{\mathit{i}}\right).\mathit{t}.{\mathit{v}}_{\mathit{y}}$$

(2)

We can now introduce the time $t_i$ taken for the action potential to reach point ${\overrightarrow{s}}_{m,i}$ from the reference ${\overrightarrow{s}}_0$ related to the position of the action potential at time 0 so that:

$$\begin{array}{c}{\mathit{x}}_{\mathit{m},\mathit{i}}={\mathit{x}}_{\mathbf{0}}+{\mathit{t}}_{\mathit{i}}\cdot {\mathbf{\upsilon }}_{\mathit{x}}\\ {\mathit{y}}_{\mathit{m},\mathit{i}}={\mathit{y}}_{\mathbf{0}}+{\mathit{t}}_{\mathit{i}}\cdot {\mathbf{\upsilon }}_{\mathit{y}}\\ {\mathit{z}}_{\mathit{m},\mathit{i}}={\mathit{z}}_{\mathbf{0}}\end{array}$$

(3)

By combining the equations above, it is possible to find the position of the minimum approach point ${\overrightarrow{s}}_{m,i}$ in terms of the velocity, the electrode coordinates (x_i,y_i) and the arbitrary reference point (x_o,y_o); the z-coordinate is invariant due to its orthogonal orientation to the plane of the electrodes.

$$\begin{array}{l}{x}_{m,i}=\frac{x_i\cdot v_x{}^2+v_x\cdot v_y(y_i-y_o)+x_o\cdot v_{y^{}}{}^2}{v_x{}^2+v_y{}^2}\\ y_{m,i}=\frac{y_i\cdot v_y{}^2+v_x\cdot v_y(x_i-x_o)+y_o\cdot v_x{}^2}{v_x{}^2+v_y{}^2}\end{array}$$

(4)

The coordinates for the electrode i and the origin are known, but not the velocity or the position of the reference point. We can use the information provided by another electrode, say electrode j; by introducing $\mathsf{\Delta }{x}_{ij}=\left(x_{m,i}-x_{m,j}\right)$, $\mathsf{\Delta }{y}_{ij}=\left(y_{m,i}-y_{m,j}\right)$ and $\mathsf{\Delta }{t}_{ij}=\left(t_i-t_j\right)$ one obtains, after some calculations:

$$\Delta t_{ij}\cdot (v_x{}^2+v_y{}^2)=v_x\cdot \Delta x_{ij}+v_y\cdot \Delta y_{ij}$$

(5)

Equation (5) has two unknown variables, v_x and v_y, and therefore one more electrode (here terms electrode k) is required. Using Equation (5) for both pairs of electrode gives the following:

$$\begin{array}{l}\frac{v_x}{v_x{}^2+v_y{}^2}=\frac{\Delta t_{ij}\cdot \Delta y_{ik}-\Delta t_{ik}\cdot \Delta y_{ij}}{\Delta y_{ij}\cdot \Delta x_{ik}-\Delta x_{ij}\cdot \Delta y_{ik}}=K_x\\ \frac{v_y}{v_x{}^2+v_y{}^2}=\frac{\Delta t_{ij}\cdot \Delta x_{ik}-\Delta t_{ik}\cdot \Delta x_{ij}}{\Delta y_{ij}\cdot \Delta x_{ik}-\Delta x_{ij}\cdot \Delta y_{ik}}=K_y\end{array}$$

(6)

It is useful at this stage to introduce the following complex number:

$$K_x-i.K_y=\frac{v_x-i.v_y}{v_x^2+v_y^2}=\frac{1}{(v_x+i.v_y)}$$

(7)

Such that:

$$v_x+i.v_y=\frac{1}{K_x-i.K_y}$$

(8)

This complex variable can provide the velocity of the action potential by taking the real and imaginary part of its inverse. The values K_x and K_y depend on two variables: the position of the electrodes, and the difference in time between the observation of the maximum signal of an AP at the two different electrodes, assuming that the probe is inserted such that the general direction of neural travel is in the direction of the x-axis, and that the difference in times of an action potential seen by two electrodes can be calculated from the recordings. It is therefore possible to process K_x and K_y and to estimate the velocity; one can also determine the position of the neuron by using the model of propagation of the electric field, since the potential due to a charge q is given by the expression

$$V_i(t)=\frac{1}{4\pi \epsilon {\epsilon }_0}\frac{q}{D_i}=\frac{q\prime }{D_i}$$

(9)

where q′ = q/4πεε₀. This expression is difficult to solve in this form, but can be simplified by using a local frame of reference {X, Y, Z} rotated by an angle $\alpha =\arctan (v_y/v_x)$, such that the neuron is aligned with the X-axis. This reduces the number of unknowns to three, the coordinates Y and Z and the charge q. This can be further simplified by considering the potential V_i for each electrode when the action potential is at the point ${\overrightarrow{s}}_{m,i}$ so that no coordinates along the X-axis appear in the equations. Under these conditions, combining Equations (1) and (9) for all the electrodes in the local frame of reference provides the following solution:

$$\{\begin{array}{l}{q}^2=\frac{\left(Y_i^2\left(Y_j^{}-Y_k^{}\right)+Y_j^2\left(Y_k^{}-Y_i^{}\right)+Y_k^2\left(Y_i^{}-Y_j^{}\right)\right)}{S_i^{-2}(Y_k^{}-Y_j^{})+S_j^{-2}(Y_i^{}-Y_k^{})+S_k^{-2}(Y_j^{}-Y_i^{})}\\ Y=\frac{1}{2}\frac{S_i^{-2}(Y_k^2-Y_j^2)+S_j^{-2}(Y_i^2-Y_k^2)+S_k^{-2}(Y_j^2-Y_i^2)}{S_i^{-2}(Y_k^{}-Y_j^{})+S_j^{-2}(Y_i^{}-Y_k^{})+S_k^{-2}(Y_j^{}-Y_i^{})}\\ Z^2=\frac{q^2}{S_i^2}-{\left(Y_{}-Y_i\right)}^2=\frac{q^2}{S_j^2}-{\left(Y_{}-Y_j\right)}^2=\frac{q^2}{S_k^2}-{\left(Y_{}-Y_k\right)}^2\end{array}$$

(10)

where the variables S_i, S_j and S_k represent the maximum amplitudes of an aAP seen by electrodes i, j and k at positions Y_i, Y_j and Y_k. The coordinates Y and Z are expressed in the local frame and must be rotated back to the initial frame.

This preliminary model is made using an assumption that the AP can be represented by a single point charge, which is a coarse estimation. It is possible to improve this by introducing the shape of the action potential. We can consider the action potential as a distribution of charges Q along the neuron, or, $Q\prime =Q/4\pi \epsilon {\epsilon }_0$ which is a function of time t and position s. Equation (9) then becomes:

$$\begin{array}{l}{V}_i(t)={\displaystyle \underset{neuron}{\mathsf{\int }}\frac{Q\prime (t-s/v)}{D_i(s)}ds}\\ ={\displaystyle \underset{neuron}{\mathsf{\int }}\frac{Q\prime (t-s/v)}{D_{m,i}\sqrt{1+{(s-s_{m,i})}^2/D_{m,i}{}^2}}ds}\end{array}$$

(11)

where $D_{m,i}$ is the minimum distance between probe i and the neuron. Mathematically, Equation (11) can be regarded as a convolution between the charge function Q′ and a function F defined by:

$$F(s)=\frac{1}{\sqrt{1+{\left(\frac{s-s_{m,i}}{D_{m,i}{}^2}\right)}^2}}$$

(12)

Then, Equation (11) becomes:

$${\mathit{V}}_{\mathit{i}}\left(\mathit{t}\right)=\frac{\mathbf{1}}{{\mathit{D}}_{\mathit{m},\mathit{i}}}\left(\mathit{Q}\mathbf{\prime }\mathit{F}\right)\left(\mathit{t}\right)=\frac{{\mathit{q}}_{\mathit{e},\mathit{i}}\mathbf{\prime }\left(\mathit{t}\right)}{{\mathit{D}}_{\mathit{m},\mathit{i}}}$$

(13)

where

$$q{\prime }_{e,i}(t)={\displaystyle \int Q\prime (t-s/v).F(s)ds}$$

(14)

Equation (13) for V_i is similar to Equation (9) for a point charge, but here the equivalent charge depends on the position of the electrode and time. This makes the analytical calculation of the electric potential much more difficult. However, using Equation (12), it is possible to avoid this problem by estimating the first-order error provided by the point-charge approximation. One can note the following:

$$\frac{\partial q_{e,i}\prime (t)}{\partial D_{m,i}}={\displaystyle \underset{s}{\mathsf{\int }}Q(t-s/v)F(s)\frac{1}{D_{m,i}}\frac{1}{1+D_{m,i}^2/{(s-s_i)}^2}ds}$$

(15)

We can set an upper bound to the value of the fraction under the integrand as follows:

$$\frac{1}{1+D_{m,i}^2/{(s-s_i)}^2}<\frac{1}{1+D_{m,i}^2/s_{\max }{}^2}$$

(16)

And then:

$$\left|\frac{\partial q{\prime }_i(t)}{\partial D_{m,i}}\right|<\frac{1}{D_{m,i}}\frac{1}{1+D_{m,i}^2/s_{\max }{}^2}\left|q{\prime }_i(t)\right|$$

(17)

s_max can be estimated from the value of the velocity and the AP duration using a plot of function F, then one can express the first-order Taylor series of the electric potential as:

$$V_i(t)\approx \frac{q\prime (t)}{D_{m,i}}\left(1+\frac{1}{D_{m,i}}\frac{1}{1+D_{m,i}^2/s_{\max }^2}\left|D_i-D\right|\right)=\frac{q\prime (t)}{D_{m,i}}\left(1+{\epsilon }_i\right)$$

(18)

The error ε_i can be processed from the estimated values in order to estimate the error in the position of the neuron.

The resolution of the system is defined by the sampling rate and electrode spacing; for example, if the signals are sampled at 40 kHz, then only differences in time higher than 25 μs can be detected; if in this instance the distance between electrodes is 80 μm then the maximum detectable velocity would be 3.2 m s⁻¹, which is sufficient for non-myelinated nervous systems, such as neural cultures or invertebrate study. A 1 MHz, 100 μm system would be required to detect the fastest 100 m s⁻¹ action potentials.

3. Methods

In order to test the model, we analyzed recordings from implantable neural probes inserted into the ventral nerve cords of desert locusts. The locust model was selected due to the fact that it allowed recordings from a large number of unmyelinated neurons organized in a structure that is both spatially limited (to the dimensions of the cord) and directionally well defined (with axons generally running along the cord axis). The approach was used with a variety of triode penetrating electrode arrays, with electrodes in triangular arrays ranging between 20 and 80 µm apart and with areas between 64 and 177 µm² [11,12,13], manufactured with recording surfaces in either gold or iridium. Past investigation of the effects of varying electrode size and positioning suggests our electrodes were suitable for positioning, being sufficiently small to avoid averaging of the peak in transit, and large enough to minimize errors due to electrode impedance [14].

The probes were mounted into a headstage amplifier with further remote amplification stages being applied before sampling; details of the recording setup can be found elsewhere [15]. The signal was acquired using an AT-MIO 16-F National Instruments acquisition (Austin, TX, USA) card at 40 kHz/channel using a LabView (National Instruments, USA) and saved for later analysis. Recordings were taken from the ventral nerve cord, metathoracic ganglia and jumping leg nerve of the desert locust (L. Migratoria). The experimental protocol for recording can be found elsewhere [15]; experiments were performed in compliance with the relevant laws and University of Surrey guidelines.

Locusts were anaesthetized, decapitated and dissected to expose the cavity where the nerve cord was found, which was filled with locust saline [12] at room temperature. The probe was inserted into the nerve using a micromanipulator such that the electrode plane was parallel to the direction of the nerve. If no spontaneous signals were observed, a gentle touch was applied to the leg, or warm saline (at 28–30 °C) was poured over the body. When the signals started appearing, recordings were taken.

The recorded signals from the locust model were input into a MatLab (the Mathworks, Natick, MA, USA) program to implement the detection routine described above. Prior to analysis, frequencies below 100 Hz were removed from the signal using a FFT function, followed by wavelet de-noising using the wavelet sym5, selected for its resemblance to an action potential. A total of eight recording sessions were performed; during each of which, significant numbers of APs were recorded. Following analysis, nerves were examined histologically and measured. An example of an AP appearing in multiple electrode channels with different amplitudes and peak times can be seen in Figure 2. Peaks were detected by thresholding, and then tracked through multiple frames, with the amplitude and duration being extracted for each electrode.

4. Results and Discussion

The implementation of the algorithm is shown in Figure 3 for a typical neural recording. On the top right the figure is a 3D diagram representing the space above the probe. The probe, as illustrated, sits on the x–y plane pointing towards the positive part of the y axis. The line crossings represent the position of the axons with color indicating the velocity of the action potentials corresponding to a given axon. The plane on top of the probe is a transverse section of the space above the probe where neurons are crossing. The box shows a section of the nerve, the top left and bottom diagrams present the different views of nerve. The top left is a superior view; the bottom left a transverse section and the bottom right a lateral view. The blue crosses represent the position of the neurons at closest approach to the probe. In the analysis, the large majority of lines are running parallel to the probe, indicating that the probe was inserted with an angle parallel to the nerve, as was indeed the case.

We found that with a small percentage (typically 5–10% across multiple experiments) were observed to lie at distances outside the estimated boundary of the ventral nerve cord as estimated using histology. We found that a similar number showed estimated velocity vectors that were at significant variance from the majority. It should be noted that in such cases, greater distance from the probe is determined by smaller signal amplitude, which in turn means the signals are more susceptible to noise.

The locust ventral nerve cord is well characterized in terms of the distributions if neurons of different diameter [15,16,17,18,19,20]. Since action potential velocity varies with neuron diameter, this gives an indication of the distribution of anticipated action potentials, and provides a method of assessing whether the system provides results in line with expectations. The distribution of measured velocities for a typical experiment can be seen in Figure 4. It should be noted that given the nature of the experiment (decapitation followed by stimulation of the leg) that we would anticipate that only afferent signals (heading to the brain) would be generated, with no efferents (in the direction of signals coming from the brain) except reflex arcs. This is in line with observations, in which 93% of signals were observed in one direction. This approach can be expanded to include AP velocity. In an average locust ventral nerve cord [19], 64% of neurons have diameters greater than 4 μm and corresponding velocities greater than 1 ms⁻¹; 24% have diameters between 2 and 3 μm and velocities of 0.5–1 ms⁻¹; and 12% have diameters below 2 μm and velocities below 0.5 ms⁻¹. Our figures showed 74% having greater than 1 ms⁻¹, 15% having velocity between 0.5 and 1 ms⁻¹; and 11% with velocities below 0.5 ms⁻¹. This is in general good agreement with expectations, but suggests the system is more sensitive to large-diameter, rapid neurons than small-diameter, slower neurons. It could be noted that visual inspection of the curve on Figure 4 suggests an unusual spike at ca. 1.2 ms⁻¹, which may be attributable to significant numbers of APs from a single neuron skewing the results, and without which the measured and estimated numbers would be more alike. Alternatively, the discrepancy may be attributed to different size ratios being used for afferents vs. efferents, since only the former will be present here (to any significant degree), which are not accounted for in the anatomical studies described above. Examination of the size distribution of a typical analysis (Figure 4b) reveals capture across a radius of approximately 100 µm, but also highlights multiple collocated spikes (indicated by darker markers). This may have an impact on the velocity distribution in Figure 4a, but also suggests the system is effective in classifying neural sources by observing multiple firings at similar special locations.

There are a number of issues that affect the accuracy of the approach. We have investigated variability according to electrode size [14] and fabrication material [21] and found that variability in these over the ranges we have considered does not significantly affect the signals acquired by these electrodes.

The proposed algorithm uses a simple model of field propagation that does not account for discrepancies that might be caused by differences in the tissue environment or the electrical properties of the extracellular medium [22,23]. As outlined above, the necessity for identifying the peak makes the system susceptible to electrical noise, particularly where the AP is relatively small (such as at a distance from the probe, or where the neural diameter is small). Whilst the effect will also be present for signals that pass near the electrode but start at a distance from it, such errors can be eliminated by identification of the signal at closest approach. Conversely, classification by monitoring the AP as it moves through space and time offers advantages over methods that relay on a “snapshot” of nearby APs; for example, whilst the shapes of APs will be different at different positions in space (electrodes). The use of “point” electrodes (ca. >20 µm across) reduces the effect of “flattening” observed with larger electrodes, making the identification of peak timing easier, whilst the taking of multiple frames of the same AP makes identification of the true “peak” detection down to a fraction of a millisecond simpler. This offers significant benefits over “traditional” spike classification studies that only use a single electrode and use amplitude and shape characteristics to sort spikes to specific neurons

However, given the reach of the algorithm (taking a snapshot across a region ca. 200 µm in diameter as seen in Figure 4b), it may be more applicable to the detection of whole nerve patterns rather than identification of individual neurons, enabling the use of pattern recognition to identify patterns relating to specific functions. Such approaches are in common use in the central nervous system for the identification of intention [24,25,26,27]; our technology would enable a similar approach to be taken in both central and peripheral nervous systems. It can also be noted that although the system presented here is designed for detection of APs in non-myelinated nervous systems, the principles of triangulation using multiple electrodes simultaneously tracking APs by amplitude and length has application in the identification of nodes of Ranvier in the region of detection for application in myelinated systems. Alternatively, the work demonstrates that multiple electrodes can be used to track multiple AP sources in, for example, a planar electrode array for in vitro neural analysis, particularly where an interlayer (such as a thin gel) could be used to increase the triangulation effect by increasing the z-dimension between electrode and neuron planes.

5. Conclusions

We have developed a model that can be used to track the movement of action potentials in a non-myelinated nervous system. The model developed has been applied to the detection of action potentials in a locust nerve, in which it has detected the positions of the neurons with some success when compared with known anatomy of the nerve. Once fully developed, such a system may also be applied to the more complex issue of tracking action potentials in myelinated nerve tissue, with evident applications in neural recording for interfacing with prostheses.