NL2036134B1 - Spike sorting under influence of external stimulation artefacts - Google Patents

Abstract

A computer-implemented method for determining a model of a plurality of neurons of a neural tissue; the method comprising: obtaining spontaneous readout data, the spontaneous readout data comprising a time series of arrays of spontaneously generated activation signals, preferably voltage signals, the spontaneously generated activation signals being generated in absence of external stimulation; obtaining evoked readout data, the evoked readout data comprising a time series of arrays of activation signals, preferably voltage signals, the activation signals being generated under influence of said external stimulation; determining an electrical imaging, El, model based on the spontaneous readout data, the El model representing activation signal propagation over the plurality of neurons in absence of said external stimulation; determining an artefact model based on the evoked readout data, the artefact model representing artefact propagation over the plurality of neurons due to said external stimulation; and determining an aggregate model based on a combination of the El model and the artefact model, the aggregate model representing activation signal propagation as well as artefact propagation over the plurality of neurons under influence of said external stimulation.

Description

SPIKE SORTING UNDER INFLUENCE OF EXTERNAL STIMULATION ARTEFACTS

TECHNICAL FIELD

The present disclosure generally relates to neural modelling. Particular embodiments relate to a computer-implemented method for determining a model of a plurality of neurons of a neural tissue, a computer-implemented method of spike sorting for a plurality of neurons of a neural tissue, a computer program, and a data processing apparatus.

BACKGROUND

Bi-directional electronic neural interfaces, which utilize electrical recording and stimulation to establish effective communication with the nervous system, have emerged as a promising avenue in the field of neurotechnology. This innovative approach enables the precise calibration of electrical stimulation while simultaneously recording the evoked neural response. This ability to facilitate bidirectional communication between electronic devices and the nervous system has far-reaching implications for various applications, including neuro-prosthetics, neural rehabilitation, and fundamental neuroscience research.

In the quest to develop reliable bi-directional electronic neural interfaces, researchers have explored several core approaches over the years. One key avenue involves the use of multi-electrode arrays (MEAs), which consist of multiple microelectrodes implanted in close proximity to neural tissue. MEAs allow for the simultaneous recording of neural activity and the delivery of electrical stimulation, thus enabling real- time feedback and adaptation of stimulation parameters.

Another established method in the field is the use of intracortical microelectrode arrays, which are particularly relevant for neural prosthetics. These arrays are designed to penetrate the cerebral cortex, providing high-resolution recordings of neural signals, and have been utilized in research aimed at developing brain-computer interfaces (BCIs) for individuals with motor impairments.

SUMMARY

Despite these advancements, the field faces challenges related to stimulation artefacts that can obscure evoked neural signals. These artefacts often arise due to the inherent nature of electrical stimulation, requiring researchers to develop sophisticated signal processing techniques and stimulation strategies to mitigate their impact. However, these techniques and strategies are cumbersome, and therefore, itis an aim of at least some embodiments according to the present disclosure to streamline the compensation of such artefacts.

Accordingly, there is provided in a first aspect of the present disclosure, a computer- implemented method according to claim 1, for determining a model of a plurality of neurons of a neural tissue; the method comprising: - obtaining, by said computer, spontaneous readout data, the spontaneous readout data comprising a time series of arrays of spontaneously generated activation signals, preferably voltage signals, the spontaneously generated activation signals being generated in absence of external stimulation; - obtaining, by said computer, evoked readout data, the evoked readout data comprising a time series of arrays of activation signals, preferably voltage signals, the activation signals being generated under influence of said external stimulation; - determining, by said computer, an electrical imaging, El, model based on the spontaneous readout data, the El model representing activation signal propagation over the plurality of neurons in absence of said external stimulation; - determining, by said computer, an artefact model based on the evoked readout data, the artefact model representing artefact propagation over the plurality of neurons due to said external stimulation; and - determining, by said computer, an aggregate model based on a combination of the

El model and the artefact model, the aggregate model representing activation signal propagation as well as artefact propagation over the plurality of neurons under influence of said external stimulation.

In other words, the method defines a computer-based analysis method that: receives spontaneous readout data (free from artefacts) and calculates an El model based on that spontaneous readout data, receives evoked readout data (incorporating artefacts) and calculates an artefact model based on that received evoked readout data, and finally combines those two models to form an aggregate model.

For a signal to be generated under influence of external stimulation may be taken to mean that the signal is generated in the vicinity and presence of such external stimulation, or that the signal is generated in the vicinity of and immediately after such external stimulation.

An artefact (defined in the feature of the “artefact propagation”) is an electrical signal that is at least partially, typically largely or fully, attributable to the external stimulation.

In an alternative but equivalent point of view, the artefact is any signal that is not consistent with one being generated by the neural substrate.

In various embodiments, the steps of determining the El model, the artefact model, and the aggregate model comprise, respectively: - determining, by said computer, an El model prescribed by a set of parameters describing the propagation of the electrical signal across electrodes while minimizing an error between collected data and an output generated by said model; - determining, by said computer, an artefact model prescribed by a set of parameters describing the propagation of the electrical signal across electrodes and another set of parameters describing an influence of said external stimulation, while minimizing an error between collected data and an output generated by said model; and - determining, by said computer, an aggregate model combining the El model and the artefact model, by adding the El model and the aggregate model so as to minimize a joint error between collected data and an output generated by said aggregate model.

In various embodiments, the step of determining the artefact model based on the evoked readout data comprises: - obtaining, by said computer, a definition of at least one array position where said external stimulation was applied to the neural tissue; and

- discarding, by said computer, any array elements from at least one time slice of the time series of the arrays, said array elements having activation signals with an amplitude exceeding a predefined sensor dynamic range threshold; and/or discarding, by said computer, at least one array element from at least one time slice of the time series of the arrays, based on the at least one array position where said external stimulation was applied.

Advantageously, by discarding the array element where the stimulation was applied, the direct impact of the stimulation on the downstream analysis procedure can be reduced, while still benefiting from measurements made at other array elements in the relative neighbourhood of the stimulation.

Said definition may e.g. be obtained from an independent entity responsible for (or at least aware of) the external stimulation, e.g. an operator controlling a stimulating electrode array, or may e.g. be obtained via signal analysis, for instance by determining that one (or a few) specific array elements are higher than some predefined threshold or are extremely high or are undefined or are defined as "Not A

Number" or something similar, and can thus be interpreted to fall outside of the dynamic range of the sensor array and therefore to represent the point(s) of the external stimulation.

In an exemplary embodiment, the array element(s) to be discarded does/do not have to remain the same permanently, but only during a specific time window. If, for instance, there are 512 array elements (because there would be 512 electrodes in the electrode array used), at the array position where said external stimulation was applied, i.e. at the stimulation site, say electrode 1, data from time 0 ms to 5 ms can be ignored, and spontaneous data can be registered at the same electrode afterwards, if the electrodes are bidirectional. Then the method may, for example, proceed similarly on electrode 2 and 3 from 2 ms to 4 ms, etc.

In various embodiments, the at least one array element to be discarded further comprises all array elements directly adjacent to the array position where said external stimulation was applied.

In a particular exemplary embodiment, array elements are considered directly adjacent to one another (in particular to the array position where said external stimulation was applied) if the distance between their centres is less than 1.5 times the pitch distance 5 of the array, preferably at most v(2) times the pitch distance.

This means that in a square grid, not only the horizontal and vertical neighbours may be included, but also the diagonal neighbours, and that in hexagonal grids all neighbours may be included.

In various embodiments, the combination is a linear combination.

This has the benefit of straightforward computation, with adequate model performance.

In various embodiments, the combination includes measurement noise.

In various embodiments, the method comprises: - arranging at least one electrode sensor array on the plurality of neurons of the neural tissue; - recording the spontaneous readout data; - providing the spontaneous readout data to said computer; - applying external stimulation to at least one neuron of the plurality of neurons of the neural tissue; -recording the evoked readout data, under influence of said applied external stimulation; and - providing the evoked readout data to said computer.

In various further developed embodiments, the at least one electrode sensor array is a bi-directional electrode sensor array.

In various further developed embodiments, the electrode sensor array is a high-density multi-electrode sensor array.

In various further developed embodiments, the electrode sensor array comprises at least 9 electrodes, preferably at least 512 electrodes, more preferably at least 65636 electrodes, at a pitch of between 10 um and 100 um, preferably between 10 um and 60 pum.

In various further developed embodiments, the step of applying the external stimulation comprises delivering a charge-balanced electrical pulse via at least one stimulation electrode of the electrode sensor array, the electrical pulse having a duration of between 50 us and 1 ms, preferably between 100 ps and 250 us.

The electrical pulse is preferably triphasic, but biphasic or any other shape may be used instead. The electrical pulse is preferably a current pulse, but may instead be a voltage pulse.

The method works, in principle, with any stimulation type, duration and waveform.

Different stimulation amplitudes (e.g. increasing over time) may preferably be used, because the threshold of each specific neuron may be unknown. However, the method is not linked to the number of stimulation amplitudes. Multiple trials are preferably used because it is a stochastic process (i.e. if you stimulate always with the same amplitude, you do not get always the same result), it is desired to build enough statistics to calculate the probability of activating the neuron, but of course the specific number of trials is not a limitation to the present invention.

In various embodiments, the plurality of neurons is part of an animal, primate or human neural tissue, preferably a retina, outside (ex vivo) or inside (in vivo) the organism.

Additionally, there is provided in a second aspect of the present disclosure, a computer-implemented method according to claim 13, of spike sorting for a plurality of neurons of a neural tissue, in particular under influence of external stimulation artefacts; the method comprising: - obtaining, by said computer, an aggregate model determined according to any preceding claim, with respect to the same plurality of neurons of the neural tissue;

- obtaining, by said computer, a definition of at least one array position where external stimulation was applied to the neural tissue; - obtaining, by said computer, newly evoked readout data, the newly evoked readout data comprising a time series of arrays of activation signals, preferably voltage signals, the activation signals being generated under influence of said external stimulation; and - estimating, by said computer, a time series of likely inputs to the plurality of neurons of the neural tissue corresponding with the newly evoked readout data, based on the aggregate model.

In other words, the aggregate model may be used to do “spike sorting”, under stimulation artefacts, i.e. to determine which inputs of firing neurons were (most) likely to have generated the observed outcome.

Additionally, there is provided in a third aspect of the present disclosure, a computer program according to claim 14, comprising instructions which, when the program is executed by a computer, cause the computer carry out the method of any one of the preceding claims.

Additionally, there is provided in a fourth aspect of the present disclosure, a data processing apparatus according to claim 15, comprising means for carrying out the method of any one of the preceding claims.

In a particular embodiment, the data processing apparatus may comprise at least one processor and at least one memory storing the above-described computer program.

The embodiments described herein are provided for illustrative purposes and should not be construed as limiting the scope of the invention. It is to be understood that the invention encompasses other embodiments and variations that are within the scope of the appended claims. The invention is not restricted to the specific configurations, arrangements, and features described herein. The invention has wide applicability and should not be limited to the specific examples provided. The embodiments disclosed are merely exemplary, and the skilled person will appreciate that various modifications and alternative designs can be made without departing from the scope of the invention.

It will be appreciated that considerations and advantages applicable to some embodiments of the invention, e.g. the method, may also be applicable to other embodiments of the invention, e.g. the computer program, mutatis mutandis, and vice versa.

BRIEF DESCRIPTION OF THE DRAWINGS

In the following description, a number of exemplary embodiments will be described in more detail, to help understanding, with reference to the appended drawings, in which:

Figure 1 schematically illustrates an exemplary embodiment of the method according to the present disclosure;

Figure 2 schematically illustrates another exemplary embodiment of the method according to the present disclosure;

Figures 3A-3D schematically illustrate another exemplary embodiment of the method according to the present disclosure;

Figures 4A and 4B schematically illustrate a further exemplary embodiment of the method according to the present disclosure;

Figures 5A and 5B schematically illustrate another exemplary embodiment of the method according to the present disclosure;

Figures 6A and 6B schematically illustrate an exemplary embodiment of the method according to the present disclosure, and show in particular an overview of the proposed spike sorting using dynamical systems and input estimator;

Figure 7 schematically illustrates an exemplary embodiment of the method according to the present disclosure, and in particular shows the performance of the input estimator when analysing synthetic data;

Figure 8 shows prediction results from the best and worst performing El models for three different electrodes;

Figure 9 schematically illustrates a SNRMSE distribution for all El models;

Figures 10A and 10B schematically illustrates prediction results of the artefact model for two electrodes at low (Figure 10A, at 0.45uA) and high (Figure 10B, at 3.71uA) stimulation amplitudes;

Figure 11 shows the ANRMSE as a function of the stimulation amplitude;

Figure 12 shows a comparison between retrieved inputs and input templates when El 1 and El 25 are fed to the input estimator;

Figure 13 schematically illustrates retrieved inputs of neuron 25 and neuron 1 when stimulation data is fed to the input estimator at low amplitude (0.4514) and high amplitude (3.71uA);

Figure 14A and 14B schematically illustrate probability of activation as a function of the stimulation amplitude for neuron 25 and neuron 1, respectively; and

Figure 15 shows the activation thresholds obtained by humans and the proposed method, i.e. a comparison of activation thresholds between human spike sorting and the proposed method.

DETAILED DESCRIPTION

In the following detailed description, specific exemplary embodiments will be described in more detail. It will be appreciated that the skilled person can transpose one or more features (and their resulting effects) described below with reference to particular embodiments having particular features to other embodiments having different features, as long as this makes technical sense to the skilled person.

Bi-directional electronic neural interfaces use electrical recording and stimulation to communicate with the nervous system, and permit accurate calibration of electrical stimulation by simultaneously recording the evoked response. However, unavoidable stimulation artefacts can obscure evoked neural signals.

The present disclosure describes a novel approach using dynamical control systems to identify electrically-evoked neural activity in the presence of electrical artefacts.

The present disclosure develops various embodiments of the method as described above, that exploit the distinct spatiotemporal propagation of neural activity and electrical artefacts to identify and distinguish spikes from distinct cells (i.e. spike sort).

This approach also makes it possible to ignore signals from electrodes severely corrupted by the stimulation artefact (typically, the stimulating electrode) while still successfully differentiating evoked spikes from electrical artefacts.

The proposed method can reliably distinguish evoked spikes from electrical artefacts by exploiting the distinct spatiotemporal propagation recorded with a dense large-scale multi-electrode array. In the future, this approach could be potentially extended to handle real-time closed-loop stimulation applications and multi-channel stimulation. 1 Introduction

Bi-directional neural interfaces play an increasingly important role in neurotechnology to communicate with the nervous system using multi-electrode arrays (MEAS). These interfaces promise to revolutionize scientific discovery and clinical therapeutics through closed-loop neuromodulation [1-5]. Specially, a BNI performs two major tasks.

On the one hand, it stimulates neurons to produce targeted patterns of neural activity that are useful for the scientific or clinical application [6-8]. On the other hand, it performs electrical recording to observe natural neural activity and to calibrate the activity evoked by the interface [9-16].

A critical challenge in recording electrically-evoked activity is that the voltage produced by injecting the current into the electrode-electrolyte impedance produces a stimulation artefact that is often large enough to obscure the evoked neural signal of interest [17]. Because of the large time constants of the electrode-tissue impedance, the artefact can last for several milliseconds after stimulation [18,19] and can thus overlap in time with evoked spikes. This substantially complicates the process of identifying and segregating spikes from different cells (spike sorting). Therefore, the artefact and the neural activity of interest must be distinguished [20, 21].

Several approaches have been proposed to use the temporal properties of spikes and artefacts to perform spike sorting [15,22-26]. In template subtraction methods, the estimated artefacts are subtracted from the measurements to isolate neural activity [19,27-29] and identify spikes [30]. However, obtaining templates of the artefact in isolation is not always possible [31].

However, relatively little has been done to exploit the distinct spatiotemporal propagation of electrical artefacts and spikes [18,32]. The present disclosure describes a novel approach using dynamical control systems to model the spatiotemporal propagation of spikes and artefacts and exploit their differences to identify evoked neural activity. Specifically, the present disclosure includes a unique dynamical model for the stimulation artefact and for each neuron recorded during spontaneous activity. Then, to identify evoked spikes after stimulation, the present disclosure includes an estimation of which combination of dynamical models (i.e. which neurons firing) was most likely to produce the recorded voltage response across all electrodes. Notably, the method does not require recordings from the stimulation electrode itself, which typically has an artefact that saturates recording electronics, enabling lower-power electronics. The effectiveness of the described approach has been tested on large scale multi-electrode ex vivo recordings from primate retina, and the results have been compared to human-supervised spike sorting.

Figure 1 schematically illustrates an exemplary embodiment of the method according to the present disclosure. The figure illustrates a bird-eye view of spike sorting in the presence of stimulation artefacts. A charge-balanced tri phasic waveform (indicated with an arrow) is used by the stimulator in this example. The neural interface records the voltage immediately after the stimulation ends. The recordings include the residual artefact(s), and evoked and/or spontaneous spikes (MEA regions of interest are as circles and polygons, respectively). Post-processing is needed to analyse the recordings and separate the effect of different phenomena (spikes and artefacts). 2 Materials and Methods 2.1 Experimental Data Description

We analysed voltage recordings from primate retinal ganglion cells (RGCs) obtained with a MEA with E = 512 electrodes at 60um pitch and a sampling frequency of 20kS/s.

To identify the characteristic signal associated with a spike in each cell, the electrical images (Els) for N = 25 neurons of interest recorded over 71 time samples (3.55ms) were examined. The El reveals the average recorded voltage trace of the spike in a particular neuron on all electrodes [14]. Typically, the voltage waveform of the spike recorded near the cell body of a given neuron Figure 2 (201) is larger than the voltage waveform of the spike propagating along its axon (202, 203).

Figure 2 schematically illustrates another exemplary embodiment of the method according to the present disclosure. The figure shows El data for neuron number 18.

On the left side of the figure, time-series are shown for three different electrodes 201, 202, 203. On the right side of the figure, an MEA heatmap 204 is shown with the neuron 205 depicted in shaded contours for the relevant electrodes. The heatmap shows the footprint of the neuron 205 on the MEA and the shading encodes the root mean square (RMS) of each electrode’s recordings when a spike from the putative neuron is present.

Electrical stimulation was performed with a brief (150us) triphasic charge-balanced current pulse delivered at the stimulating electrode (39 stimulus amplitudes, 25 trials for each amplitude). The voltage after the stimulation ended was recorded in 55 samples on all electrodes (2.75ms). Typically, the stimulating electrode (Figure 3A, 3B) recorded an electrical artefact much larger than non-stimulating electrodes (Figure 3C, 3D). The artefact duration was approximately 0.4ms on non-stimulating electrodes, but nearly 2ms on the stimulating electrode, comparable to or greater than the latency of directly evoked spikes (0.4ms - 0.6ms) [18]. The magnitude of the artefact increased with stimulus amplitude on all electrodes, but the stimulating electrode exhibited discontinuities in the artefact amplitude range due to the stimulation hardware design [17].

Figures 3A-3D schematically illustrate another exemplary embodiment of the method according to the present disclosure. The figures show recordings of a stimulating electrode 301 and two non-stimulating electrodes 302, 303for different stimuli amplitudes:

Figure 3A shows the position 300 of the electrodes on the MEA.

Figure 3B corresponds with the stimulating electrode 301.

Figure 3C corresponds with the non-stimulating electrode 302.

Figure 3D corresponds with the non-stimulating electrode 303.

Note that the range of amplitude for the stimulating electrode 301 is of course much greater than the range of amplitude for the non-stimulating electrodes 302, 303, and note also that the relatively closer non-stimulating electrode 302 shows a higher amplitude than the relatively further non-stimulating electrode 303, in relation to the array position of the stimulating electrode 301. 2.2 Dynamical Systems Approach

We use dynamical systems to model the spatiotemporal evolution of the voltage recorded by the MEA in the presence of spikes from one or more neurons and in the presence of a stimulation artefact [33, 34]. 2.2.1 El Models

The El model for a given neuron n is a set of equations (e.g. eq. 1) that describes the propagation of the voltage across all E = 512 electrodes when the neuron fires a spike.

The model consists of a predefined input vector u, € R2that initiates the dynamical system, the state of the neuron x’, € RE" at each point in time defined over a subset of the electrodes relevant for the cell (E, < E, see region 205 in Figure 2), and two matrices A, and B that describe the evolution of x!, over time. Finally, the output yn €

RE indicates the voltage contribution of neuron n to all electrodes on the array (the output is set to zero for the electrodes that are not relevant for the cell). Thus: x= Apxin + Bath, y= Cun (1)

The crucial matrix A, defines how the state of the neuron at the next time step depends on the current state, and captures the spatiotemporal correlation between the electrodes. The matrix B, defines how the next state depends on the input, and C, defines how the output depends on the current state.

The model parameters are learned such that if the input templates r, are injected, the model generates outputs (i.e. the voltage on each electrode) that approximately match the El of neuron n. The input templates are two unit-magnitude pulses that trigger rising and falling phases of the spike, respectively.

Figures 4A and 4B schematically illustrate a further exemplary embodiment of the method according to the present disclosure. The figure shows input templates and their effect on the El model. In Figure 4A, the circles show the effect of the input templates on the electrodes, as defined by matrix B. The arrows show the spatial propagation of the model across electrodes, as defined by matrix A,. Figure 4B shows that the input template includes two pulses 401 and 402. 2.2.2 Artefact Model

The artefact mode! describes the propagation of voltage across the electrodes caused by artefacts due to electrical stimulation. Similarly to the El model, the output ya € RE indicates the voltage recorded on all electrodes immediately after stimulation. The empirically observed artefact propagates radially outward from the stimulating electrode and decays over space. Hence, only the closest electrodes E, < E to the stimulating electrode will record the artefact and are considered as the states of the model, namely x's € RE? at time t. Moreover, the stimulating electrode shows discontinuities with respect to the stimulus amplitude due to the design of the stimulation system (see Figure 3B). Thus, we discard the stimulating electrode from the state variables. As with the El model, y/ais equal to state variable x', for electrodes

Es, and zero for the other electrodes. As with the El model, the artefact model is learned such that an approximation to the artefact is observed in the output ya if the input u's has the template "4 € R3 (Figure 5B) linearly scaled by the stimulus amplitude g. Hence, for the stimulus amplitude gq, the input to the model is u's = gra. In summary, the dynamical model of the artefact is given by the following equation, referenced as equation (2): x= Aaxla + Balla,

Vaz Caxa where the vector u's € R%is the input of the model, the Ba matrix defines how the next state depends on the input of the model, the matrix As defines how the next state depends on the current state, and the matrix Ca defines how the output depends on the current state. A,, Ba and C, are obtained from the average of stimulation data across 25 trials.

Figures 5A and 5B schematically illustrate another exemplary embodiment of the method according to the present disclosure. The figures show input templates and their effect on the artefact model, i.e., the figure demonstrates how the artefact model characterizes the propagation of the artefact on MEAs. Figure 5A shows dynamics of the artefact. In Figure 5A, the circles show the effect of the template inputs on the electrodes, as defined by matrix B. The arrows show the spatial propagation of the model across electrodes, as defined by matrix A,. Figure 5B shows input templates r, of artefact and shows that the input template includes three pulses 501, 502 and 503. 2.2.3 Aggregate Model

The aggregate model describes the propagation of the voltage across electrodes in the presence of both stimulation artefact and neurons firing. This model is a linear combination of the El model for all neurons, the artefact model and the measurement noise (superposition assumption in MEA recordings [23]). As shown in Figure 6A, the aggregate model outputs the electrode measurements for multiple input templates.

Specifically, the observed outputs on the electrodes are caused by injecting templates of the different sub models (El and/or artefact models) in the aggregate model. Thus, the aggregate model is given by:

MARL B iL BU, y- CX + W (3) where x= [x7y, ..., Xv, X75 ]7 is the aggregate state i = [u7y, . . . , Un ]7is the aggregate input of the El models, 4 is the input of the artefact model, and w' € Rf denotes the measurement noise [A B, B' C] are defined such that the aggregate output y satisfies y <PVper Vin + Vat W (Db) Input estimator scheme.

Figures 6A and 6B schematically illustrate an exemplary embodiment of the method according to the present disclosure, and show in particular an overview of the proposed spike sorting using dynamical systems and input estimator. In Figure 8A, it is assumed that the measured output is the superposition of the output of dynamical models for different Els and artefacts as driven by predefined input templates. The combination of the El models and artefact models is the aggregate model. In

Figure 6B, a method is illustrated to estimate the set of input templates that led to a specific measurement (see details in Section 2.2.4). The input estimator requires the electrode measurements immediately after stimulation, the aggregate model, and the list of corrupted electrodes to ignore (typically, the stimulating and nearby electrodes shown darker in the grid). Since we know a priori that stimulation happened, the artefact template is automatically estimated. By inferring which inputs are active (i.e. they are estimated to be equal to the input template), we can identify which neuron(s) are firing; hence, we can perform spike sorting. 2.2.4 Input Estimator

The aggregate model described in Section 2.2.3 determines how the electrodes’ voltage is generated based on the input sequence (Figure 6A). Conversely, the input estimator leverages the aggregate model to infer (or, equivalently, estimate) the input sequence that likely generated the observed data (Figure 6B). Notably, the estimated input sequence can be compared to the El templates to perform spike detection and sorting - i.e. if u, = r,, then the neuron n is firing a spike.

The input estimator utilizes the output data y**t=[(A7, ..., (y"*) 717 within a forward window L to estimate 4’ of i, and xX’ of x. Also, it uses the artefact inputs u's as known inputs since the timing of the stimulation is known. Additionally, / indicates the indices for the corrupted electrodes that we desire to exclude from the input estimator (i.e. the stimulating electrode in our case). Therefore, the input estimation procedure can be described as (a x" = fy, 1 ua, / (4) where the function Ay“, #71, u's, I) is the so-called input estimator.

To illustrate and validate the performance of the input estimator for spike sorting, we first consider an example using synthetic data for a single neuron firing, generated by feeding the template for the El model of one cell to the aggregate model. The generated data is then passed to the input estimator to retrieve which neuron(s) fired.

As expected, the estimated input template completely overlaps with the input template of the single El model used to generate the synthetic data (Figure 7). Hence, the input estimator can infer from the synthetic data that the neuron was firing a spike. Section 3.3 validates the input estimator with real noisy data and Section 3.4 validates the spike sorting method based on the input estimator with real data contaminated by the stimulation artefact.

Figure 7 schematically illustrates an exemplary embodiment of the method according to the present disclosure, and in particular shows the performance of the input estimator when analysing synthetic data. The input templates of El 1 are injected into the aggregate model to generate the synthetic data. The input estimator retrieves inputs that completely overlap with the input templates used to generate the synthetic data (shown on the right). 3 Results

This section shows the performance of El and artefact models, and the input estimator.

First, we provide evidence that the El models are suitable for representing real data accurately. Next, we show the performance of the artefact model for different stimulus amplitudes. We then provide evidence of the input estimator’s ability to predict the correct input sequence with spontaneous activity (i.e. without stimulation artefact) and evoked activity (i.e. with stimulation artefact). Finally, we demonstrate the effectiveness of the proposed spike sorting for different stimulation amplitudes and compare the results to human-supervised spike sorting. 3.1 El Models

To analyse the performance of the El models, we can compare the output of the El model (i.e. the predicted voltage for all electrodes) with the real measurements of the

El for all neurons. To quantify the quality of the proposed El models, we can use spike normalized root mean square error (SNRMSE): the root mean square of the difference between the real output (‘True’) and estimated output (‘Predicted’) divided by the root mean square of the El signal.

Figure 8 shows prediction results from the best and worst performing El models for three different electrodes 801, 802, 803, and 811, 812, 813. The spatial map of the El is shown on the MEA heatmap. The prediction results are for the neurons with the lowest (Figure 8B, showing prediction results for El number 9) and the highest (Figure 8A, showing prediction results for El number 11) SNRMSE. The real/true output is shown with references 821, 831, 841, 851, 861, and the estimated/predicted output is shown with references 822, 832, 842, 852, 862.

The generated output from the El models follows the El original data for the different electrodes. Notice that the electrode near the soma shows a lower SNRMSE than the other electrodes along the axon. The higher error for the neuron in Figure 8B may be caused by the location of the corresponding neuron at the edge of the MEA where there are fewer electrodes to capture its dynamical evolution. The right-skewed

SNRMSE distribution across all neurons (Figure 9) provides evidence that most El models incur SNRMSE less than the average.

Figure 9 schematically illustrates a SNRMSE distribution for all EI models. The MEA heatmap shows the neurons’ location based on their SNRMSE. The average SNRMSE is depicted with a dashed line. 3.2 Artefact Model

To analyse the performance of the artefact model, we compare the model output to the real measurements for the electrodes close to the stimulating electrode (except for the stimulating electrode itself, which is discarded from our model) for different stimulus amplitudes. Figure 10 shows how the output generated by the simulated model follows the real measurement for two electrodes. At low stimulation amplitude (Figure 10A), the model follows the average stimulation data for the electrodes. At high stimulation amplitude (Figure 10B), the model exhibits higher error in a specific interval, possibly due to the presence of an evoked spike superimposed in the recording. This possibility is consistent with the fact that higher stimulation amplitude usually leads to a higher probability of neuron activation.

Figures 10A and 10B schematically illustrates prediction results of the artefact model for two electrodes at low (Figure 10A, at 0.45uA) and high (Figure 10B, at 3.71uA) stimulation amplitudes. The circle indicates the position of the stimulating electrode.

The model output is compared with the average stimulation data across all 25 trials.

For the high amplitude cases, a putative elicited spike is shown in a dashed line.

To quantify the performance of the model, we can use artefact normalized root mean square error (AN RMSE): the root mean square of the difference between the real output and estimated output divided by the root mean square of the artefact signal.

Figure 11 shows the ANRMSE as a function of the stimulation amplitude. The

ANRMSE is highest for low stimulation amplitudes, when the artefact is small, and the output is dominated by noise in Figure 11. However, ANRMSE is reduced for higher stimulation amplitudes since the artefact becomes the dominant contributor to the output, showing that the model is capable of predicting the artefact shape. 3.3 Input Estimator

The input estimator can retrieve the input sequences based on the electrodes’ measurements. By comparing the input sequence to the input templates for each neuron, we can retrieve which neuron fired a spike, i.e. perform spike sorting. In

Section 2.2.4, we showed that the input estimator can retrieve the correct input sequence from synthetic data generated by the aggregate model (Figure 7). Here, we show that the input estimator can still recover the correct input sequence from real El data (Figure 12, showing a comparison between retrieved inputs and input templates when El 1 and El 25 are fed to the input estimator). In particular, when analysing data from EI 1 (i.e. data where only neuron 1 is firing a spike), the input estimator retrieves the input for neuron 1 with a similar shape to its predefined template, suggesting that neuron 1 is firing a spike. In contrast, the estimated input for El 25 is zero, indicating that neuron 25 is not firing a spike. Additionally, when analysing data from El 25 (i.e. data where only neuron 1 is firing a spike), the input estimator retrieves only the input template from neuron 25 and zero from neuron 1. Notice that this analysis is for measurements without a stimulation artefact.

Now, let us consider measurements immediately after stimulation that contain a stimulation artefact and can contain one or more spikes (Figure 13). In this example, the stimulation electrode is near neuron 25 and farther from neuron 1. Hence, it is expected that electrical stimulation is more likely to elicit a spike in neuron 25 than in neuron 1. At low stimulation amplitude, neither neuron fires a spike. At high stimulation amplitudes, however, the retrieved input sequence for neuron 25 suggests that it fired a spike, while still no spike was retrieved for neuron 1.

Figure 13 schematically illustrates retrieved inputs of neuron 25 and neuron 1 when stimulation data is fed to the input estimator at low amplitude (0.4514) and high amplitude (3.71uA). Input templates for each neuron are shown in dashed lines.

Neuron 25 is located behind the stimulating electrode (indicated by a circle), and neuron 1 is further away. Using a similarity function, a spike is detected only for neuron 25 at high stimulation amplitude (bottom-right panel). 3.4 Spike Sorting with Stimulation Artefact

To perform spike sorting, the input sequence estimated by the input estimator is compared to the input template for each neuron using a similarity function. If its value exceeds a predefined threshold, then we infer that a spike occurred for that specific neuron. By extending this procedure to all stimulation amplitudes over multiple trials, we obtain the activation curve for each neuron-electrode pair: the probability of a neuron firing as a function of the stimulus amplitude applied on a given electrode.

Figures 14A and 14B show an example of the activation curve of neurons 25 and 1 when stimulating using an electrode close to neuron 25 for both the proposed and human spike sorting. The proposed method results in a similar activation threshold to the human result for neuron 25 (defined as the amplitude for which the spiking probability is 50%) and no activation for neuron 1. The threshold is extracted by fitting a sigmoid curve to the scatter plot. Although the proposed method finds approximately the same threshold as human spike sorting, it presents certain differences in the sigmoid fit: (1) non-zero probability for low amplitudes, (2) higher variability, and (3) shallower curvature. In the future, adding a prior an monotonicity during spike sorting could improve these aspects.

Figure 14A and 14B schematically illustrate the probability of activation as a function of the stimulation amplitude for neuron 25 and neuron 1, respectively. The stimulation threshold (shown with a dashed line) is defined as the amplitude for which the probability of activation is 50% using a sigmoid fit of the scatter plot. The spatial map of the neurons are shown in the MEA, and a blue circle indicates the position of the stimulating electrode.

To show an overview of the performance of the proposed spike sorting, we analysed 10 different datasets in which there are 5 neurons stimulated by 10 different stimulating electrodes. Figure 15 shows the activation thresholds obtained by humans and the proposed method, i.e. a comparison of activation thresholds between human spike sorting and the proposed method. The dashed line shows the identity line. The stimulation thresholds estimated using the approach closely tracked those obtained with manual analysis (R295.1%, 10 neuron electrode pairs). 4 Conclusion 4.1 Spatio-temporal Spike Sorting

This present disclosure proposes to exploit the different spatiotemporal characteristics of spikes and electrical stimulation artefact (Figures 4a, 5a) to identify electrically evoked spikes in neural recordings. The present disclosure contrasts with prior art approaches that only exploit the different latency and amplitude of spikes and artefacts, failing to leverage the distinct spatiotemporal progression captured by high- density large-scale MEAs [19, 25]. We provided evidence that these spatiotemporal characteristics can be effectively modelled using the dynamical systems in (1) and (2), and that an input estimator can be used to retrieve the correct input sequence from

MEA recordings. The activation thresholds obtained for 10 electrode-neuron pairs were similar to those obtained with human spike sorting (Figure 15).

4.1.1 Exploiting redundancy in high-density recordings

The proposed method can also exploit the redundancy in high-density large-scale MEA recordings to remove data from contaminated electrodes while still recovering the spatiotemporal characteristics needed to perform spike sorting. The method made it possible to ignare the stimulating electrode because it presents the largest artefact and instead exploits information available on other electrodes. This approach can increase the power and area efficiency of future hardware implementations since the analog front-end does not need to record the large artefact. Furthermore, our approach overcomes the need for modelling the artefact in the stimulating electrode, which is a non-trivial task and depends on the stimulation setup. For instance, [18] needs a different model for the stimulating electrode because it presents discontinuities and is dramatically different than the other electrodes.

Since the proposed dynamical model has more outputs than inputs, the inputs could be retrieved even in the absence of one or more electrodes. Hence, this approach could generalize to different scenarios where more than one electrode is contaminated. 4.2 Caveats 4.2.1 Model Requirements and Limitations

The proposed approach relies on predefined input templates that are used to initiate the dynamical system and generate an output that resembles the El of a given neuron,

Figure 4b. Additionally, the artefact model requires input templates that may fit a specific hardware system. These input templates can however be flexibly designed and adapted to match the characteristics of a specific hardware system.

We have shown that the El models can faithfully represent the El recordings with small errors (Figure 9). However, for the neurons located on the edges of the MEA, the model may be less accurate. This is likely due to incomplete information on the El spatiotemporal propagation, which generates discontinuities in the MEA recordings.

4.2.2 Computational Efficiency

The low computational complexity for spike sorting is important in real-time and closed-loop applications [35]. For instance, the work proposed in [18, 23, 36-38] uses computationally complex matching pursuit methods for spike sorting, where whole recordings are compared to the spiking templates of each cell. In contrast, our approach only exploits a short time window of recordings to recover input templates via the input estimator. Hence, the input estimator functions as a computationally efficient filter over the recordings similar to the Wiener filters in [27, 31]. 4.2.3 Background Noise

Our approach relies on the availability of Els for each neuron obtained by the spike sorting method in [23] in the absence of stimulation. The present disclosure does not focus on low signal-to-noise ratio (SNR) data because Els with a peak amplitude less than 30uV are discarded. Thus, the discarded spikes are present in the background noise and may worsen the performance of the approach [39]. At the same time, the background noise is coloured with spatiotemporal dependencies [40, 41], and the methods in [23, 42] can estimate its covariance. Therefore, it provides a chance to research the correlation in the background noise that may be useful in reducing the range of SNR that can be accommodated by the algorithm.

Another issue for the proposed algorithm is bundle activation, which is the result of stimulating a bundle of axons belonging to cells whose somas are outside the region covered by the MEA [43]. These spikes are also part of the background noise and may worsen the performance of the proposed method. Furthermore, if the goal is to precisely stimulate neurons at single cell and cell type resolution, then bundle activation needs to be avoided. This can be done by combining the proposed method with the method developed in [44] that detects bundle activation based on its spatiotemporal propagation characteristics. 4.2.4 Artefacts

In modelling the artefact, a crucial concern is that the artefact signal is unsupervised and unknown beforehand. The stimulation data consists of the evoked neural activity, the artefact, and the background noise. Thus, we have to recover artefacts from the stimulation data to fit the artefact model. To remove the effect of noise and spontaneous neural activity, we used the average of the stimulation data from different trials similar to [18]. This works well for low-amplitude stimulation because it does not systematically lead to evoked spikes. However, as shown in Figure 14A, stimulation with high amplitudes leads to a higher probability of neuron activation. Hence, the artefact cannot be estimated by averaging for high stimulation amplitude. In literature, to model the artefact, many approaches rely on assumptions on artefact timing, lack of saturation, linearity, and decline with distance from the stimulating electrode [32, 45, 46]. Here, we simplify this task by proposing a method that does not need to model the artefact in the stimulating electrode where it is most significant and can be nonlinear. 4.2.5 Hyperparameters

The models, the input estimator, and the spike detector all require hyperparameters.

For the models, we experimentally fitted the hyperparameters to minimize the

SNMRSE and AMRSE. For the input estimator, we tuned the hyperparameters based on the performance when estimating synthetic and real data (Figure 7 and Figure 12).

For the spike detector, we relied on human annotations for adjusting hyperparameters such as detection thresholds. Hence, we require multiple human-annotated labels, which are costly and time-consuming [18, 47]. 4.3 Towards closed-loop multi-channel stimulation

We conjecture that spatiotemporal spike sorting is suitable to enable real-time stimulation and recording at single-cell resolution in closed-loop applications. Towards this goal, future work should focus on automating the proposed method (e.g. for the hyper-parameter tuning) and designing a real-time hardware implementation.

In addition, we believe that the proposed framework can cope with multi-channel stimulation scenarios. This could have a large impact in basic neuroscience and clinical applications since multi-channel stimulation enables an effective strategy to spatially control the spiking of multiple neurons [48, 49]. To address the multi-channel stimulation scenario, the proposed method could factor in multiple artefact models in the aggregate model and retrieve the input templates of Els in the presence of different artefacts simultaneously.

As used in this application and in the claims, the singular forms “a,” “an,” and “the” include the plural forms unless the context clearly dictates otherwise. The systems, apparatus, and methods described herein should not be construed as limiting in any way. Instead, the present disclosure is directed toward all novel and non-obvious features and aspects of the various disclosed embodiments, alone and in various combinations and sub-combinations with one another. The disclosed systems, methods, and apparatus are not limited to any specific aspect or feature or combinations thereof, nor do the disclosed systems, methods, and apparatus require that any one or more specific advantages be present or problems be solved. Any theories of operation are to facilitate explanation, but the disclosed systems, methods, and apparatus are not limited to such theories of operation.

Although the operations of some of the disclosed methods are described in a particular, sequential order for convenient presentation, it should be understood that this manner of description encompasses rearrangement, unless a particular ordering is required by specific language set forth below. For example, operations described sequentially may in some cases be rearranged or performed concurrently. Moreover, for the sake of simplicity, the attached figures may not show the various ways in which the disclosed systems, methods, and apparatus can be used in conjunction with other systems, methods, and apparatus. Additionally, the description sometimes uses terms like “obtaining” and “outputting” to describe the disclosed methods. These terms are high-level abstractions of the actual operations that are performed. The actual operations that correspond to these terms will vary depending on the particular implementation and are readily discernible by the skilled person.

It will be appreciated that numerous specific details are set forth in order to provide a thorough understanding of the examples described herein. However, it will be understood by the skilled person that the examples described herein can be practiced without these specific details. In other instances, methods, procedures and components have not been described in detail so as not to obscure the related relevant feature being described. The drawings are not necessarily to scale and the proportions of certain parts may be exaggerated to better illustrate details and features. The description is not to be considered as limiting the scope of the examples described herein.

Reference throughout this specification to “one example,” “an example,” “certain examples,” or "exemplary implementation” means that a particular feature, structure, or characteristic described in connection with the feature and/or example may be included in at least one feature and/or example of claimed subject matter. Thus, the appearances of the phrase “in one example,” “an example,” “in certain examples,” “in certain implementations,”, “in various embodiments”, or other like phrases in various places throughout this specification are not necessarily all referring to the same feature, example, and/or limitation. Furthermore, the particular features, structures, or characteristics may be combined in one or more examples and/or features.

References

[1]D. G. Muratore and E. J. Chichilnisky, Artificial Retina: A Future Cellular-Resolution

Brain Machine Interface. Cham: Springer International Publishing, 2020, pp. 443-465.

[2] M. W. Slutzky, “Brain-machine interfaces: Powerful tools for clinical treatment and neuroscientific investigations,” The Neuroscientist, vol. 25, pp. 139-154, 4 2019.

[3] M. A. Lebedev and M. A. Nicolelis, “Brain-machine interfaces: past, present and future,” TRENDS in Neurosciences, vol. 29, no. 9, pp. 536-546, 2006.

[4] D. Lam, H. A. Enright, J. Cadena, S. K. G. Peters, A. P. Sales, J. J. Osburn, D. A.

Soscia, K. S. Kulp, E. K. Wheeler, and N. O. Fischer, “Tissue-specific extracellular matrix accelerates the formation of neural networks and communities in a neuron-glia co-culture on a multi-electrode array,” Scientific Reports, vol. 9, p. 4159, 12 2019.

[5] M. David-Pur, L. Bareket-Keren, G. Beit-Yaakov, D. Raz-Prag, and Y. Hanein, “All- carbon nanotube flexible multi-electrode array for neuronal recording and stimulation,”

Biomedical Microdevices, vol. 16, pp. 43-53, 2 2014.

[6] L. H. Jepson, P. Hottowy, K. Mathieson, D. E. Gunning, W. Dabrowski, A. M. Litke, and E. J. Chichilnisky, “Focal electrical stimulation of major ganglion cell types in the primate retina for the design of visual prostheses,” Journal of Neuroscience, vol. 33, pp. 7194-7205, 4 2013.

[7] G. K. Moghadam, R. Wilke, G. J. Suaning, N. H. Lovell, and S. Dokos, “Quasi- monopolar stimulation: A novel electrode design configuration for performance optimization of a retinal neuroprosthesis,” PLoS ONE, vol. 8, 8 2013.

[8] G. A. Goetz and D. V. Palanker, “Electronic approaches to restoration of sight,”

Reports on Progress in Physics, vol. 79, p. 096701, 9 2016.

[9] K. A. Moxon and G. Foffani, “Brain-Machine Interfaces beyond Neuroprosthetics,”

Neuron, vol. 86, no. 1, pp. 55-67, 2015.

[10] M. A. Lebedev and M. A. L. Nicolelis, “Brain-machine interfaces: From basic science to neuroprostheses and neurorehabilitation,” Physiological Reviews, vol. 97, pp. 767-837, 4 2017.

[11] Y. Yang, S. Qiao, O. G. Sani, J. |. Sedillo, B. Ferrentino, B. Pesaran, and M. M.

Shanechi, “Modelling and prediction of the dynamic responses of large-scale brain networks during direct electrical stimulation,” Nature Biomedical Engineering, vol. 5, pp. 324-345, 4 2021.

[12] R. Segev, J. Goodhouse, J. Puchalla, and M. J. Berry, “Recording spikes from a large fraction of the ganglion cells in a retinal patch,” Nature Neuroscience, vol. 7, pp. 1154-1161, 10 2004.

[13] G. Portelli, J. M. Barrett, G. Hilgen, T. Masquelier, A. Maccione, S. D. Marco, L.

Berdondini, P. Kornprobst, and E. Sernagor, “Rank order coding: A retinal information decoding strategy revealed by large-scale multi-electrode array retinal recordings,” eNeuro, vol. 3, pp. 844-853, 2016.

[14] A. M. Litke, N. Bezayiff, E. J. Chichilnisky, W. Cunningham, W. Dabrowski, A. A.

Grillo, M. Grivich, P. Grybos, P. Hottowy, S. Kachiguine, R. S. Kalmar, K. Mathieson,

D. Petrusca, M. Rahman, and A. Sher, “What does the eye tell the brain?:

Development of a system for the large-scale recording of retinal output activity,” J/EEE

Transactions on Nuclear Science, vol. 51, pp. 1434-1440, 8 2004.

[15] F. Franke, M. Natora, C. Boucsein, M. H. Munk, and K. Obermayer, “An online spike detection and spike classification algorithm capable of instantaneous resolution of overlapping spikes,” Journal of Computational Neuroscience, vol. 29, pp. 127-148, 8 2010.

[16] S. M. Potter, A. E. Hady, and E. E. Fetz, “Closed-loop neuroscience and neuroengineering,” Frontiers in Neural Circuits, vol. 8, 2014.

[17] P. Hottowy, A. Skoczen, D. E. Gunning, S. Kachiguine, K. Mathieson, A. Sher, P.

Wiacek, A. M. Litke, and W. Dabrowski, “Properties and application of a multichannel integrated circuit for low artefact, patterned electrical stimulation of neural tissue,”

Journal of Neural Engineering, vol. 8, no. 6, p. 066005, Nov. 2012.

[18] G. E. Mena, L. E. Grosberg, S. Madugula, P. Hottowy, A. Litke, J. Cunningham,

E. J. Chichilnisky, and L. Paninski, “Electrical stimulus artefact cancellation and neural spike detection on large multi electrode arrays,” PLOS Computational Biology, vol. 13, no. 11, p. 1005842, Nov. 2017.

[19] T. Hashimoto, C. M. Elder, and J. L. Vitek, “A template subtraction method for stimulus artefact removal in high-frequency deep brain stimulation,” Journal of

Neuroscience Methods, vol. 113, pp. 181-186, 2002.

[20] M. Pais-Vieira, A. P. Yadav, D. Moreira, D. Guggenmos, A. Santos, M. Lebedev, and M. A. Nicolelis, “A closed loop brain-machine interface for epilepsy control using dorsal column electrical stimulation,” Scientific Reports, vol. 6, 9 2016.

[21] J. Müller, D. Bakkum, and A. Hierlemann, “Sub-millisecond closed-loop feedback stimulation between arbitrary sets of individual neurons,” Frontiers in Neural Circuits, 12 2012.

[22] D. J. Caldwell, J. A. Cronin, R. P. Rao, K. L. Collins, K. E. Weaver, A. L. Ko, J. G.

Ojemann, J. N. Kutz, and B. W. Brunton, “Signal recovery from stimulation artefacts in intracranial recordings with dictionary learning,” Journal of Neural Engineering, vol. 17, 4 2020.

[23] J. W. Pillow, J. Shlens, E. J. Chichilnisky, and E. P. Simoncelli, “A model-based spike sorting algorithm for removing correlation artefacts in multi-neuron recordings,”

PLOS ONE, vol. 8, no. 5, pp. 1-14, 05 2013.

[24] J. Lee, C. Mitelut, H. Shokri, I. Kinsella, N. Dethe, S. Wu, K. Li, E. B. Reyes, D.

Turcu, E. Batty, Y. J. Kim, N. Brackbill, A. Kling, G. Goetz, E. Chichilnisky, D. Carlson, and L. Paninski, “Yass: Yet another spike sorter applied to large-scale multi-electrode array recordings in primate retina,” bioRxiv, 2020. [Online]. Available: https://www.biorxiv.org/content/early/2020/03/20/2020.03.18.997924

[25] C. Sekirnjak, P. Hottowy, A. Sher, W. Dabrowski, A. M. Litke, and E. J.

Chichilnisky, “Electrical stimulation of mammalian retinal ganglion cells with multi- electrode arrays,” Journal of Neurophysiology, vol. 95, pp. 3311-3327, 6 20086.

[26] ——, “High-resolution electrical stimulation of primate retina for epiretinal implant design,” Journal of Neuroscience, vol. 28, pp. 4446-4456, 4 2008.

[27] M. Schelles, J. Wouters, B. Asamoah, M. M. Laughlin, and A. Bertrand, “Objective evaluation of stimulation artefact removal techniques in the context of neural spike sorting,” Journal of Neural Engineering, vol. 19, 2 2022.

[28] A. E. Mendrela, J. Cho, J. A. Fredenburg, V. Nagaraj, T. I. Netoff, M. P. Flynn, and E. Yoon, “A bidirectional neural interface circuit with active stimulation artefact cancellation and cross-channel common-mode noise suppression,” IEEE Journal of

Solid-State Circuits, vol. 51, pp. 955-965, 4 2016.

[29] S. Culaclii, B. Kim, Y.-K. Lo, L. Li, and W. Liu, “Online artefact cancelation in same-electrode neural stimulation and recording using a combined hardware and software architecture,” IEEE Transactions on Biomedical Circuits and Systems, vol. 12, no. 3, pp. 601-613, 2018.

[30] M. S. Lewicki, “A review of methods for spike sorting: The detection and classification of neural action potentials,” Network: Computation in Neural Systems, vol. 9, 1998.

[31] M. S. Najafabadi, L. Chen, K. Dutta, A. Norris, B. Feng, J. W. Schnupp, N.

Rosskothen-Kuhl, H. L. Read, and M. A. Escabi, “Optimal multichannel artefact prediction and removal for neural stimulation and brain machine interfaces,” Frontiers in Neuroscience, vol. 14, 7 2020.

[32] D. J. O'Shea and K. V. Shenoy, ‘“Eraasr: An algorithm for removing electrical stimulation artefacts from multi-electrode array recordings,” Journal of Neural

Engineering, vol. 15, 2 2018.

[33] O. Nelles, Nonlinear System Identification. Springer Berlin Heidelberg, 2001.

[34] R. Isermann and M. Münchhof, Identification of Dynamic Systems. Springer Berlin

Heidelberg, 2011.

[35] N. P. Shah and E. J. Chichilnisky, “Computational challenges and opportunities for a bi-directional artificial retina,” Journal of Neural Engineering, vol. 17, no. 5, Oct. 2020, publisher: IOP Publishing Ltd.

[36] S. Garcia, A. P. Buccino, and P. Yger, “How do spike collisions affect spike sorting performance?” Neuroscience, preprint, Dec. 2021.

[37] M. Pachitariu, N. A. Steinmetz, S. N. Kadir, M. Carandini, and K. D. Harris, “Fast and accurate spike sorting of high-channel count probes with kilosort,” in Advances in

Neural Information Processing Systems, D. Lee, M. Sugiyama, U. Luxburg, |. Guyon, and R. Garnett, Eds., vol. 29. Curran Associates, Inc., 2016. [Online]. Available: https://proceedings.neurips.cc/paper_files/paper/2016/file/1145a30ff80745b56fb0cec f65305017-Paper.pdf

[38] P. Yger, G. L. Spampinato, E. Esposito, B. Lefebvre, S. Deny, C. Gardella, M.

Stimberg, F. Jetter, G. Zeck, S. Picaud, J. Duebel, and O. Marre, “A spike sorting toolbox for up to thousands of electrodes validated with ground truth recordings in vitro and in vivo,” eLife, vol. 7, p. e34518, Mar. 2018.

[39] M. Sahani, “Latent Variable Models for Neural Data Analysis,” Ph.D. dissertation,

California Institute of Technology, 1999.

[40] M. S. Fee, P. P. Mitra, and D. Kleinfeld, “Automatic sorting of multiple unit neuronal signals in the presence of anisotropic and non-gaussian variability,” Journal of Neuroscience Methods, vol. 69, no. 2, pp. 175-188, 1996.

[41] Z. Yang, Q. Zhao, E. Keefer, and W. Liu, “Noise characterization, modeling, and reduction for in vivo neural recording,” in Advances in Neural Information Processing

Systems, Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, Eds., vol. 22. Curran Associates, Inc., 2009, p. 2180-2168.

[42] F. Franke, R. Quian Quiroga, A. Hierlemann, and K. Obermayer, “Bayes optimal template matching for spike sorting — combining fisher discriminant analysis with optimal filtering,” J Comput Neurosci, vol. 38, no. 3, pp. 439-459, Jun. 2015.

[43] L. E. Grosberg, K. Ganesan, G. A. Goetz, S. S. Madugula, N. Bhaskhar, V. Fan,

P. Li, P. Hottowy, W. Dabrowski, A. Sher, A. M. Litke, S. Mitra, and E. J. Chichilnisky, “Activation of ganglion cells and axon bundles using epiretinal electrical stimulation,”

Journal of Neurophysiology, vol. 118, no. 3, pp. 1457-1471, Sep. 2017.

[44] P. Tandon, N. Bhaskhar, N. Shah, S. Madugula, L. Grosberg, V. H. Fan, P.

Hottowy, A. Sher, A. M. Litke, E. J. Chichilnisky, and S. Mitra, “Automatic Identification of Axon Bundle Activation for Epiretinal Prosthesis,” /EEE Transactions on Neural

Systems and Rehabilitation Engineering, vol. 29, pp. 2496-2502, 2021, publisher:

Institute of Electrical and Electronics Engineers Inc.

[45] A. Zhou, B. C. Johnson, and R. Muller, “Toward true closed-loop neuromodulation: artefact-free recording during stimulation,” Current Opinion in Neurobiology, vol. 50, pp. 119-127, 6 2018.

[46] J. M. Weiss, S. N. Flesher, R. Franklin, J. L. Collinger, and R. A. Gaunt, “Artefact- free recordings in human bidirectional brain-computer interfaces,” Journal of Neural

Engineering, vol. 16, no. 1, Feb. 2019, publisher: Institute of Physics Publishing.

[47] A. H. Barnett, J. F. Magland, and L. F. Greengard, “Validation of neural spike sorting algorithms without ground-truth information,” Journal of Neuroscience

Methods, vol. 264, pp. 65-77, 2016. [Online]. Available: https://www.sciencedirect. com/science/article/pii/S0185027016300036

[48] S. Lee, J. Park, J. Kwon, D. H. Kim, and C.-H. Im, “Multi-channel transorbital electrical stimulation for effective stimulation of posterior retina,” Sci Rep, vol. 11, no. 1, p. 8745, May 2021.

[49] L. H. Jepson, P. Hottowy, K. Mathieson, D. E. Gunning, W. Dabrowski, A. M. Litke, and E. J. Chichilnisky, “Spatially Patterned Electrical Stimulation to Enhance

Resolution of Retinal Prostheses,” J. Neurosci., vol. 34, no. 14, pp. 4871-4881, Apr. 2014.

Claims (15)

CONCLUSIONS 1. A computer-implemented method for determining a model of a plurality of neurons of a neural tissue; the method comprising: - obtaining, by the computer, spontaneous readout data, the spontaneous readout data comprising a time series of arrays of spontaneously generated activation signals, preferably voltage signals, the spontaneously generated activation signals being generated in the absence of external stimulation; - obtaining, by the computer, evoked readout data, the evoked readout data comprising a time series of arrays of activation signals, preferably voltage signals, the activation signals being generated under the influence of the external stimulation; - determining, by the computer, an electrical imaging model, E1, based on the spontaneous readout data, the E1 model representing propagation of the activation signal over the plurality of neurons in the absence of the external stimulation; - computing an artifact model based on the elicited readout data, the artifact model representing artifact propagation over the plurality of neurons due to the external stimulation; and - computing an aggregated model based on a combination of the E1 model and the artifact model, the aggregated model representing propagation of the activation signal as well as artifact propagation over the plurality of neurons under the influence of the external stimulation. 2. A method according to claim 1, wherein the steps of determining the El model, the artifact model, and the aggregated model comprise, respectively: - determining, by the computer, an El model prescribed by a set of parameters describing the propagation of the electrical signal over electrodes while minimizing an error between collected data and an output generated by said model; - determining, by the computer, an artifact model prescribed by a set of parameters describing the propagation of the electrical signal over electrodes and another set of parameters describing an influence of the external stimulation, while minimizing an error between collected data and an output generated by said model; and - determining, by the computer, an aggregated model as a combination of the El model and the artifact model, by bringing together the El model and the aggregated model so as to minimize a joint error between collected data and an output generated by said aggregation model. 3. A method according to any preceding claim, wherein the step of determining the artifact model based on the evoked readout data comprises: - obtaining, by the computer, a definition of at least one array position where the external stimulation was applied to the neural tissue, and - omitting, by the computer, any array elements from at least one time portion of the time series of the arrays, which array elements represent activation signals with an amplitude exceeding a predefined dynamic range threshold of the sensor; and/or omitting, by the computer, at least one array element from at least one time portion of the time series of the arrays, based on the at least one array position where the external stimulation was applied. 4. A method according to any preceding claim, wherein the at least one array element intended to be omitted further comprises all immediately adjacent array elements to the array position at which the external stimulation was applied. 5. A method according to any preceding claim, wherein the combination is a linear combination. 6. A method according to any preceding claim, wherein the combination contains measurement noise. 7. A method according to any preceding claim, comprising: - applying an electrode sensor array to the plurality of neurons of the neural tissue; - making a recording of the spontaneous readout data; - providing the spontaneous readout data to the computer; - applying external stimulation to at least one neuron of the plurality of neurons of the neural tissue; - making a recording of the evoked readout data, under the influence of the applied external stimulation; and - providing the evoked readout data to the computer. 8. The method of claim 7, wherein the at least one electrode sensor array is a bidirectional electrode sensor array. 9. The method of any one of claims 7 to 8, wherein the electrode sensor array is a high-density multi-electrode sensor array. 10. Method according to any one of claims 7 to 9, wherein the electrode sensor array comprises at least 9 electrodes, preferably at least 512 electrodes, more preferably at least 85636 electrodes, with a spacing of between 10 um and 100 um, preferably between 10 um and 60 um. A method according to any one of claims 7 to 10, wherein the step of applying the external stimulation comprises delivering a charge-balanced electrical pulse via at least one stimulation electrode of the electrode sensor array, said electrical pulse having a duration of between 50 µs and 1 ms, preferably between 100 ms and 250 µs. 12. A method according to any preceding claim, wherein the plurality of neurons is part of an animal, primate or human nervous tissue, preferably a retina, outside (ex vivo) or inside (in vivo) the organism. 13. A computer-implemented method for peak sorting for a plurality of neurons of a neural tissue, the method comprising: - obtaining, by the computer, an aggregated model determined according to any one of the preceding claims, with respect to the same plurality of neurons of the neural tissue; - obtaining, by the computer, a definition of at least one array position where external stimulation was applied to the neural tissue; - obtaining, by the computer, newly evoked readout data, the newly evoked readout data comprising a time series of arrays of activation signals, preferably voltage signals, the activation signals being generated under the influence of the external stimulation; and - estimating, by the computer, a time series of probable inputs to the plurality of neurons of the neural tissue corresponding to the newly evoked readout data, based on the aggregated model. 14. Computer program comprising instructions which, when executed by a computer, cause the computer to perform the method according to any one of the preceding claims. 15. Data processing device comprising means for carrying out the method according to any one of claims 1 to 13.