Ricardo Escolá

Christophe Pouzat

Antoine Chaffiol

Blaise Yvert

Isabelle E Magnin

Régis Guillemaud

SIMONE: a realistic neural network simulator to reproduce MEA-based recordings

SIMONE: A Realistic Neural Network Simulator to

Reproduce MEA-Based Recordings Ricardo Escolá, Christophe Pouzat, Antoine Chaffiol, Blaise Yvert, Isabelle E. Magnin, and Régis Guillemaud Abstract-Contemporary multielectrode arrays (MEAs) used to record extracellular activity from neural tissues can deliver data at rates on the order of 100 Mbps. Such rates require efficient data compression and/or preprocessing algorithms implemented on an application specific integrated circuit (ASIC) close to the MEA. We present SIMONE (Statistical sIMulation Of Neuronal networks Engine), a versatile simulation tool whose parameters can be either fixed or defined by a probability distribution. We validated our tool by simulating data recorded from the first olfactory relay of an insect. Different key aspects make this tool suitable for testing the robustness and accuracy of neural signal processing algorithms (such as the detection, alignment, and classification of spikes). For instance, most of the parameters can be defined by a probabilistic distribution, then tens of simulations may be obtained from the same scenario. This is especially useful when validating the robustness of the processing algorithm. Moreover, the number of active cells and the exact firing activity of each one of them is perfectly known, which provides an easy way to test accuracy.

Synaptic weight between postsynaptic neuron and presynaptic neuron .

Correction factor for amplitude spatial attenuation.

I. INTRODUCTION

T HE strain put on multielectrode array (MEA) data acqui- sition systems has recently increased with the capacity to manufacture MEAs containing several hundreds of electrodes [START_REF] Olsson | A three-dimensional neural recording microsystem with implantable data compression circuitry[END_REF], [START_REF] Moulin | Contribution à létude et à la réalisation dun système électronique de mesure et excitation de tissu nerveux à matrices de microélectrodes[END_REF]. Larger arrays should moreover become available in the coming years. This current and anticipated augmentation in the number of recording sites implies an increase in the incoming amount of data. As each electrode delivers approximately 15 KB per second, a 1000-electrode MEA would require a 120 Mbps1 data path to empty the acquisition buffer. This is clearly a problem since biomedical-compatible wireless transmission does not exceed 2 Mbps [START_REF] Wise | Wireless implantable microsystems: High-density electronic interfaces to the nervous system[END_REF]. One solution to this problem is the development of an embedded preprocessing stage able to efficiently compress data and/or extract useful information from it. An obvious issue when developing such preprocessing algorithms is the choice of the set of signals used to validate them. To this end it is customary to rely on different experimental data sets that depict "many" representative states of a neural network. However, such an approach presents an important drawback: one can never be sure of the exact spiking activity of every single neuron in the tissue ("ground truth data" [START_REF] Smith | A tool for synthesizing spike trains with realistic interference[END_REF]). Alternative schemes such as accurate simulations using Hodgkin-Huxley equations are highly resource-demanding when the number of simulated neurons is large. On the other hand, different simulation tools already exist (GENESIS [START_REF] Bower | The Book of GENESIS: Exploring Realistic Neural Models with the GEneral NEural SImulation System[END_REF] and NEURON [START_REF] Hines | The neuron simulation environment[END_REF], [START_REF] Hines | NEURON: A tool for neuroscientists[END_REF] being the most widely used). Those tools have three main limitations: 1) network architectures are hand made (large network simulations are hard to setup); 2) physiological models are deterministic (only one outcome may be expected for any chosen set of parameters); 3) complicated biological lingo is often used (making it difficult for nonexperts to fully understand how to interact with these tools).

Simpler network-level simulations using integrate-and-fire models [START_REF] Delorme | Spikenet: An event-driven simulation package for modelling large networks of spiking neurons[END_REF], [START_REF] Rochel | An event-driven framework for the simulation of networks of spiking neurons[END_REF] are computationally less demanding but do not model the most preeminent feature of neuronal activity viewed by extracellular electrodes: the spike waveforms. Furthermore, remixing experimental data with noise [START_REF] Obeid | Evaluation of spike-detection algorithms fora brain-machine interface application[END_REF]- [START_REF] Rutishauser | Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo[END_REF] is a somewhat limited approach because it only covers a reduced set of simulations, which closely mimick the extracellular recordings used as reference. Other solutions involving statistical firing rates are also possible. A very interesting approach is proposed by Smith and Mtetwa [START_REF] Smith | A tool for synthesizing spike trains with realistic interference[END_REF] in which the noise and the neural activity are correlated. In the near future we will broaden our work towards the comprehension of network dynamics. For this matter, statistical spike trains are not fully suitable.

To address these issues we propose SIMONE2 (Statistical sIMulation Of Neuronal networks Engine): a simulation engine combining low computational cost with the ability to simulate spike waveforms [START_REF] Escolá | SI-MONE: A new simulation engine of realistic MEA-based extracellular measurements[END_REF]. We worked with the environment [START_REF] Stett | A Language and Environment for Statistical Computing R Development Core Team[END_REF], an open source version of the well-known language [START_REF] Chambers | Programming with Data[END_REF]. With this engine users can: 1) simulate a given network of interconnected neurons according to a fixed or random topology, 2) compute the subthreshold intracellular potential of each neuron with the classical integrate-and-fire equations while spike waveforms are generated from templates, and 3) compute the corresponding extracellular potential recorded by a MEA.

This paper is divided into three parts. Section II, which includes a description of the physiological and mathematical models behind SIMONE, along with the experimental methods performed to collect the in vivo recordings. In Section III, we compare real and simulated data and provide an example of a full spike processing chain. Section IV presents the conclusions.

II. METHODS

A. Model Description

The simulation scenario describes a MEA and a neural network. Parameters are entered via the command-line interface or via a script file. The user can define different distributions, functions and parameters to describe the desired neural network and the associated MEA. The overall model may be fully parametrized using stochastic functions.

Three acquisition functions are available to accurately depict the electrode behavior: a time-dependent gain , which may reproduce small variations in electrode acquisition properties (such as ageing, artifacts, and gain noise); the range of sensitivity of the tip (typically 100-140 m range [START_REF] Obeid | Evaluation of spike-detection algorithms fora brain-machine interface application[END_REF]); and the time-dependent extracellular noise. These functions may evolve during simulation runtime. Automatic positioning of the electrodes is performed using three parameters: the number of electrodes per dimension (2-D or 3-D), the distance between neighboring electrodes and the coordinates of the first site.

A biological neural network is automatically generated from user-defined parameters whereby neurons are positioned in a plane and characterized along with their synaptic connections. Cells are distributed in a very thin volume. In this paper, we consider an equivalent tissue plane [see Fig. 7(d)]. A 2-D MEA is placed perpendicular to the cells. A 3-D matrix may also be created by the user. Preliminary tools are also provided to simulate 3-D configurations (network MEA).

Our dynamic model3 consists of two coupled stages (see Fig. 1): 1) intracellular simulation, based on the integrate-and-fire model and 2) extracellular simulation, based on a modified version of the classic current monopole model.

1) Physiological Model: Neural network features must be quantified in terms of three basic biological parameters: size, synaptic propagation speed, and characteristic conductivity. Moreover, several statistical distributions may also be entered by the user (otherwise, standard distributions are used). These laws define the spatial distribution of neurons; the probability for a neuron to be excitatory (inhibitory otherwise); the probability of a synapse being established between any pair of neurons (this is generally dependent on the distance between neurons and their type) and the synaptic weight; the statistical distribution of the resistance and the capacitance in the membrane model; and the threshold, inverse and resting voltages in the spiking model. Five functions control the dynamics of each neuron. First, there is a kernel to modulate the synaptic weight decay (efficiency) between two connected neurons (i.e., the instantaneous pulse and time course of postsynaptic currents due to presynaptic spikes). The time course of a postsynaptic current represents the probability of neurotransmitter release [START_REF] Uteshev | A vision of synaptic transmission between central neurons[END_REF]. A second kernel modulates spike shapes according to the interspike interval [START_REF] Quirk | Interaction between spike waveform classification and temporal sequence detection[END_REF], [START_REF] Rutishauser | Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo[END_REF]. A third function models as the time-dependent self-interaction filter of the neuron [START_REF] Paninski | Maximum likelihood estimation of a stochastic integrate-and-fire neural encoding model[END_REF]. is a feedback current that increases the amplitude of the intracellular potential thus the threshold is surpassed more easily (drastically reducing the interspike interval during a limited period of time provoking a burst). Other time-dependent functions reproduce synaptic noise and allow the user to introduce an artificial stimulation current into the model . A gallery of intracellular action potential templates must be fully designed by the user. If the number of templates is smaller than the total number of neurons, new templates are automatically created using a random linear combination of a randomlychosen pair of spikes. Templates are stored in separate files. Length, shape and sampling frequency of spike models are selected by the user.

Intracellular simulation is based on the leaky integrate-and-fire model with active conductance [START_REF] Gerstner | Spiking Neurons Models[END_REF] [see Fig. 1(a)]. In the subthreshold domain, the intracellular potential for a particular neuron is characterized via an ordinary differential (1), describing the membrane current as the sum of leakage and input currents [START_REF] Olsson | A three-dimensional neural recording microsystem with implantable data compression circuitry[END_REF] We use Euler's method to solve this equation. This is a safe choice because the membrane time constant in the subthreshold regime is typically ten to a hundred times longer than the quantization step. We put a strong emphasis on modeling the input currents: presynaptic (excitatory and inhibitory), artificial (presynaptic stimulation), membrane , noise , and self-induced (originally described by Paninski et al. [START_REF] Paninski | Maximum likelihood estimation of a stochastic integrate-and-fire neural encoding model[END_REF]). This approach allows us to depict coherent behaviors and gain deeper insight into neuronal dynamics.

Leakage current is associated with

Feedback current is related to the spiking activity of the target neuron and it may be used to produce bursting activity and other firing regimes [START_REF] Wise | Wireless implantable microsystems: High-density electronic interfaces to the nervous system[END_REF] is an artificial stimulation current defined by the user. It may be used to reproduce a "real" stimulation of the neuron, as well as spontaneous activity (thus ensuring the continuous activity of the network).

is a user-defined synaptic noise current used to overcome the lack of distant presynaptic neurons (not effectively simulated) which also contributes (weakly) to neuron dynamics. Positive (or negative) values represent excitatory (or inhibitory) neurons. The synaptic current is related to presynaptic spiking activity. The synaptic kernel can be interpreted as "the time course of a postsynaptic potential evoked by the firing of a presynaptic neuron" [START_REF] Gerstner | Spiking Neurons Models[END_REF]. Synaptic conductivity is governed by the weighted sum of presynpatic firing rates. Note that represents the decay shape of each individual synaptic current, while fixes the correct units and magnitudes (4)

(5)

When the intracellular potential reaches the threshold voltage an intracellular spike is inserted into the waveform. For action potentials to be realistic as well, we use a unique intracellular spike template assigned to each neuron 4 , attenuated according to previous firing history [START_REF] Harris | Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements[END_REF] in order to determine the spike to be used at each threshold crossing [see [START_REF] Hines | The neuron simulation environment[END_REF]]. Real action potentials can be easily embedded into the model while providing a low-cost computational solution. This is much simpler than a complete Hodgkin and Huxley model [START_REF] Hines | The neuron simulation environment[END_REF] The output of this stage is the evolution of the intracellular potential for each neuron.

2) Extracellular Acquisition Model: Extracellular potentials measured by each electrode are computed from intracellular voltages. Although a detailed approximation to the source of the action potential may involve a quadrupolar approximation [START_REF] Plonsey | Bioelectric Phenomena[END_REF], for simplicity we use a monopole model [see Fig. 1(b)]. More accurate estimations could also be made by using an extracellular model based on a 1-D linear source approximation [START_REF] Gold | On the origin of the extracellular action potential waveform: A modeling study[END_REF]. In order to be able to test the influence of different types of spatial attenuation, we introduce a factor . The extracellular voltage is inversely proportional to the distance between the neuron and the theoretical electrode raised to the power and directly proportional to the membrane current. For any given electrode , the contributions of each neuron are summed at its tip. Using a slightly modified equation for the current monopole, we end up with [START_REF] Hines | NEURON: A tool for neuroscientists[END_REF]. For the ideal monopole model equals zero. We assume each neuron is a small spherical sink5 of (membrane) current within a homogeneous environment (the neural tissue) with conductivity . It is, therefore, represented by a point source of magnitude in the simulator [START_REF] Hines | NEURON: A tool for neuroscientists[END_REF] The acquisition model for each electrode is assumed to be identical and quasi-ideal (they are "ideal" in the sense of their frequency response: all pass). We know that micro electrodes input impedance is strongly capacitive. However, forms a series circuit with the acquisition system input impedance . Since is several orders of magnitude larger than within the frequency range of interest (up to 10 kHz), may be neglected. Moreover, the equivalent gain (electrode-acquisition system couple) is constant within this frequency range. Modeling of the extracellular noise is straightforward: the user just needs to input the characteristics of the noise that have actually been measured. Finally, the simulation step must be selected. A simple rule of thumb is to select a frequency value close to the sampling frequency used in the actual experiment.

B. Remarks Concerning the Parametrization Process 1) Parameters:

It is clear that a detailed analysis of the effect and sensitivity of each parameter on the overall model would largely surpass the length of the present work. However it is important to note that a simulation mainly revolves around four important axes, each of which depends on a reduced set of parameters.

1) Individual firing rate such as spontaneous spiking and bursts are governed by membrane characteristics, and nonsynaptic input currents , that is, those currents that are not induced by presynaptic spiking activity.

2) Network dynamics and network synchronization are governed by the parameters related to synaptic connectivity ( , and ) as well as to those parameters regulating the number of neurons (density and tissue dimensions).

3) Extracellular spike waveforms perceived by the different electrodes are influenced by attenuation due to tissue conductance-related properties and dynamic modulation of spike waveforms . This is particularly important when testing spike sorting techniques. 4) Signal-to-noise ratio (SNR) is determined by explicit noise properties and dynamic electrode gain characteristics . For a more detailed explanation of the effect of each parameter, see Table I.

2) How to Test Spike Processing Algorithms: By varying physiological parameters, one may test network dynamics and how different processing algorithms respond to certain regimes and/or to distinct types of cell. Moreover, the user may also try different electrode configurations and characteristics. One possible test could be performed by studying how the algorithm performance evolves with electrode ageing or acquisition noise. Site separation may also be optimized by studying the accuracy of a tetrode-based spike sorting algorithm (or more generally a -tipped polytrode). In general, the SNR may be modified at will to easily test the number of false positives and false negatives due to a given spike detection strategy. Spike shapes may be distorted to evaluate how well a certain method achieves source separation. This is all possible because we control the scenario parameters and the simulation output (for instance the specific spike "signature" and dynamics for each neuron as well as noise levels).

C. Experiments

As proof of concept, we simulate a common experiment involving multichannel silicon microprobe arrays (NeuroNexus Technologies, Ann Arbor, MI) and an insect antennal lobe [START_REF] Christensen | Multi-unit recordings reveal context-dependent modulation of synchrony in odor-specific neural ensembles[END_REF] using SIMONE. Our insect of study is the American cockroach.

1) American Cockroach Olfactory System: Olfaction plays a major role in insect reproductive success and survival. Insects possess sensitive chemosensory systems that can detect and discriminate among a huge pool of chemicals. Odors are detected by olfactory receptor neurons (ORNs) on the antenna. The axons of ORNs project into the antennal lobe of the brain where they form synaptic contacts with other neurons in the specialized neuropilar region, the glomeruli. The invertebrate antennal lobe serves as the first stage of the central olfactory processing pathway, and is similar to the vertebrate olfactory bulb as both possess glomeruli as their characteristic functional subunits. Within the insect glomerular neuropil, the afferent receptor axons form synapses with both local interneurons (LNs, inhibitory) and projection neurons (PNs, excitatory). Then PNs convey information to higher brain centers in the mushroom bodies and the lateral protocerebrum.

The American cockroach (Periplaneta americana) is one of the few insects extensively studied in the context of olfactory processing. This is mainly due to its extraordinary endurance to electrophysiological experiments and the easy accessibility to the antennal lobe.

The antennal lobe of Periplaneta americana is innervated by receptor neurons from its associated antenna [START_REF] Ernst | A neuroanatomical study on the organization of the central antennal pathways in insects[END_REF]. Around 125 eggshaped glomeruli (100 40 m) [START_REF] Boeckh | Synaptic organization and development of the antennal lobe in insects[END_REF] are the site of contact between the ORNs, and around 700 LNs and 250 PNs of the antennal lobe for each deutocerebrum [START_REF] Ernst | A neuroanatomical study on the organization of the central antennal pathways in insects[END_REF], [START_REF] Malun | Connections between the deutocerebrum and the protocerebrum, and neuroanatomy of several classes of deutocerebral projection neurons in the brain of male Periplaneta americana[END_REF].

When the insect antenna is presented with a chemical stimulus, each responsive PN produces a temporally complex phasic response that generally outlasts the stimulus [START_REF] Laurent | Temporal representations of odors in an olfactory network[END_REF], [START_REF] Lemon | Responses of cockroach antennal lobe projection neurons to pulsatile olfactory stimuli[END_REF].

The spike trains recorded from the antennal lobe PNs carry sensory information to higher brain centers, and are of particular importance to understand the first processing step in olfaction.

2) In Vivo Recordings: Only Periplaneta americana adult males are used as experimental animals. Cockroaches are reared in an incubator with free access to food and water, at temperature around 25 C. Cockroaches are cold-anesthetized prior to the experiment. Wings, legs, and some mouth parts are removed. Each insect is restrained in an acrylic glass holder, with its head fixed with dental wax (see Fig. 2). Antennae are protected with plastic tubes around them. A window of head cuticle is opened, the tracheae on the anterior face of the brain and the sheath surrounding the antennal lobes are removed. The esophagus is cut to reduce brain movement. Fresh cockroach saline is superfused on the brain.

Recordings are made in the antennal lobe using a probe made of two silicon shanks with 16 recording electrodes 6 regrouped in four tetrodes. The shanks width is around 80 m, the separation between the midline of the shanks is 150 m, the diagonal between two consecutive sites of the same tetrode is 25 m, and the distance between two tetrodes is 150 m. The probe is gently inserted into the antennal lobe such that the two tetrodes at the tip of the probe are roughly 100 m below the surface. Signals are sampled at 12.8 kHz, amplified and band-pass filtered between kHz using an IDAC2000 amplifier and the Autospike 2000 acquisition program. 7This MEA approach allows us to record many neurons from each electrode, which implies that the collected data reflect a mixture of spikes coming from many neurons, so spike sorting analysis is performed later using SpikeOMatic [START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF], [START_REF] Pouzat | Improved spikesorting by modeling firing statistics and burst-dependent spike amplitude attenuation: A markov chain monte carlo approach[END_REF]. Moreover, each neuron can be seen by several electrodes which can improve the efficiency of spike sorting techniques.

In the Periplaneta americana antennal lobe, preliminary results (Chaffiol-Pouzat, unpublished data) showed that only PNs generate "large" action potentials, therefore LNs activity is neglected when sorting spikes in our data.

D. Model Parametrization Strategy

1) Neural Network Parameters: Recorded tissue samples are assumed to be 250m-sided square surfaces for a set of four shanks (one tetrode). Tissue is considered to be homogeneous and purely resistive (delay-free). To estimate antennal lobe conductivity, we fit the transfer function [START_REF] Hines | NEURON: A tool for neuroscientists[END_REF] to attenuate spikes of approximately 250 V at a distance of 30 m between the cell and the electrode, and to 40 V at the limit of the range of sensitivity (100 m). For that, a conductivity of 0.04 , with a heterogeneity correction of 0.5 is needed [see Fig. 5(d)]. This outcome matches nicely with the fact that conductivity of white matter in the human brain has already been reported to be around 0.05 S/m [START_REF]Radiation protection standard-Exposure limits for electric & magnetic Fields-0 hz to 3 khz Australian Radiation Protection Nuclear Safety Agency[END_REF].

American cockroach antennal lobe cells may be in contact with each other. Since they have an approximate diameter of 20 m, we can consider that the equivalent spherical sources that represent neurons are separated from each other (minimal interneuron distance) by a distance of one diameter ( m). As stated in Section II-C-1, the contribution of local (inhibitory) neurons are several times smaller than those of projection (excitatory) cells, therefore we may assume the American cockroach antennal lobe to be formed exclusively by the excitatory neurons. We also consider the neurons to be homogeneously distributed. Since each electrode measures the activity of 4-7 neurons within a 100m radius 8 , it is straightforward to calculate that neuron (PNs) density is roughly neurons/ m . We assume random synaptic probability in the American cockroach antennal lobe. That is, for every pair of ordered (presynaptic and postsynaptic) neurons, the probability of finding a synapse is 50%. In addition, synaptic speed is taken to be equal and constant for all synapses: 50 m/ms.

Concerning membrane characteristics, we were not able to find specific data regarding PNs from the antennal lobe of Periplaneta americana. However, Christensen et al. [START_REF] Christensen | Ionic currents in neurones cultured from embryonic cockroach brains[END_REF] have studied 10-15m-soma diameter neurones cultured from embryonic cockroach brains. Patch clamp measurements revealed a typical membrane resistance ( ) and time constant (

) of 850 M and 20 ms, respectively. In the absence of further data, we assume that the membrane time constant and specific resistance (i.e., resistance per unit area M m ) values will be similar for both types of neurons (since they are independent of the cells' geometry). Taking a mean cell somata radius of 6. in these three voltages are negligible, and these variations do not affect the overall dynamics of the network, we shall assume them to be equal and constant for all neurons.

Synaptic kernel is generally modeled after an exponential decay with time constant [START_REF] Gerstner | Spiking Neurons Models[END_REF]. We make an estimation of using in vivo recordings from embryonic central neurons in Drosophila, where the effect of presynaptic spikes gradually vanishes at approximately 15 ms [START_REF] Baines | Electrophysiological development of central neurons in the drosophila embryo[END_REF]. We can safely estimate this decay with being 5 ms (one third of the elapsed time), see Fig. 3. We normalize this function to unity and make the synaptic weight carry the actual magnitude of the postsynaptic effect [START_REF] Obeid | Evaluation of spike-detection algorithms fora brain-machine interface application[END_REF] The synaptic weight quantifies the efficiency of a synapse. As no precise data is available concerning this parameter for Periplaneta americana, we make a first estimation based on work published by Stricker and Redman describing evoked responses of CA1 pyramidal rat cells [START_REF] Stricker | Quantal analysis based on density estimation[END_REF]. For rat cells, they report that mean peak excitatory postsynaptic current (EPSC) is around 5 . EPSC may be expressed as the synaptic current response for a single presynaptic action potential, when . We can now estimate the mean synaptic weight to be given by (11) [START_REF] Rutishauser | Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo[END_REF] Experiments show that the modulation of the spike waveform does not depend on the previous spiking history [START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF]. However, a dispersion in the amplitude is inherent to acquisition. Measurements show that amplitude dispersion from the typical spike template is symmetrically distributed % around the corresponds to the typical spontaneous spiking frequency for each neuron. Spike waveforms are not shown for simplicity. mean action potential amplitude. This variation does not depend on interspike intervals. Nevertheless, we can "trick" SI-MONE to modulate spikes using the bursting kernel even though this function is -dependent. We define as for every interspike interval. We could choose to include an -related component, but we rather prefer not to do so, as bursts are rarely noticed in our recordings. Moreover, since this is a proof of concept demonstration, we restrict our simulations to simple behaviors (no bursting, no firing rate adaptation, etc.). For the same reason, the self-interaction current is neglected: equals 0. The stimulation current is used to keep the tissue active, even when there is no synaptic input. Experiments show that neurons within the American cockroach antennal lobe spike spontaneously at frequencies between 4 and 20 Hz. We shall assume that 1) are uniformly distributed among neurons, between 4-20 Hz; 2) is constant for each neuron but in practice it presents a little physiological dispersion. To reproduce this behavior we excite the integrate-and-fire circuit with a Dirac pulses sequence, with impulses having a "mean" appearance frequency of . Each pulse makes the neuron fire once immediately. This causes the neuron to spike roughly spikes per second, without disturbing its stationary response. In Fig. 4, we show a typical simulated spontaneous spiking activity, induced by , for neurons 1-5 during a 1-s time window.

We assume that synaptic noise current is due to more distant loosely-correlated neurons. We chose it to be normally distributed around 20% of critical input current (the minimal current needed to spike) [START_REF] Escolá | SI-MONE: A new simulation engine of realistic MEA-based extracellular measurements[END_REF] with a standard deviation of 15 pA. Synaptic noise has the same distribution for every neuron.

The intracellular spike templates, used in the simulator, are obtained using twelve real extracellular sorted spikes 9 from the mentioned experiment. Each extracellular spike is up sampled, integrated over time and rescaled to have a maximum at 30 mV and a starting point at mV (corresponding to intracellular threshold overpassing) (see Fig. 5). A summary of neural network parameters is displayed in Table I.

coming

2) Extracellular Acquisition Parameters: Electrodes 10 are assumed to be displayed in a parallel-to-the-tissue plane. We consider a four-tipped Michigan tetrode, where electrodes are arranged at the corners of a 106m-sided square. The square is centered on top of the neural network and separated by a 20m gap of extracellular tissue [Fig. 7(c)]. Electrodes present a 100m range of sensitivity (ros). Moreover, the acquisition system for every electrode provides unity gain ( equals 1).

Experiments show [START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF] that extracellular recordings are autocorrelated at approximately the same frequency as that of spike waveforms (around 1 kHz). This is due to the synchronization of small groups of neurons. Overall noise and signal (the events) are statistically independent and linearly summed.

We extract these characteristics using an autoregressive model which is fully derived from real data. First, we calculate the autocorrelation of recorded data [see dashed curve in Fig. 6(a)]. Yule-Walker equations [START_REF] Kay | Modern Spectral Estimation[END_REF] allow us to calculate the autoregression coefficients and residuals . We use the Akaike Information Criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF] to chose the order (40th) of the autoregressive model. The residuals distribution (i.e., the error) is proved to be approximately Gaussian with 0 V mean 10 All electrodes are assumed to be identical. is mostly the noncorrelated portion of the raw signal, we model it with the random variable [see [START_REF] Stett | A Language and Environment for Statistical Computing R Development Core Team[END_REF]]. This extraction method is valid even if the signal of interest has few or no spikes at all. In such cases, error presents approximately the same distribution as the raw signal [START_REF] Stett | A Language and Environment for Statistical Computing R Development Core Team[END_REF] Finally, the simulation rate is fixed at 12.8 kHz to mimic that of the real experiment. A summary of extracellular acquisition parameters is displayed in Table I.

E. Methods Used To Analyze the Data

We compare two sets of data: 1) real and 2) simulated [electrode 1 in Fig. 7(c)] without any a priori knowledge. Unbiased (without previous data processing) metrics are used for both datasets: 1) mean and standard deviation, 2) SNR, 3) quantiles, 4) autocorrelation and power spectral density. The raw signal is formed by spikes and noise. As stated before, they are independent. SNR is the power ratio between a signal (spikes) and the background (extracellular) noise [see [START_REF] Chambers | Programming with Data[END_REF]]. Noise may be extracted for each data series using the residuals of the autoregressive model (see Section II-D-2). We estimate signal power as the square of the root mean square (rms) value for all data points above three times the standard deviation of the noise. As the noise has a Gaussian-shaped distribution, its power is the square of the standard deviation of the residuals for each dataset [START_REF] Chambers | Programming with Data[END_REF] where is the power.

Finally, we use SpikeOMatic [START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF] to illustrate how spikes can be effectively classified. We treat each simulated recording site individually (Table II). First, we apply a prefiltering stage to every recording site. For simplicity, we use a low-pass box filter of length three. The spike detection threshold is fixed at five times the median absolute deviation of the whole recording. Both local maxima (positive peaks) and local minima (negative valleys) are detected. A 2.5-ms-window is cut around single extremum events. Spikes are aligned in such a way that absolute maximas (peaks or valleys) are positioned at 1 ms within the extraction window. Principal component analysis [START_REF] Jolliffe | Principal Component Analysis[END_REF] is applied to reduce spike dimensionality (2.5 ms 12.8 kHz means 32 points) to four principal components. Spikes are finally classified using K-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] to identify the neuron that generated each event (spike sorting). Because each electrode records the activity of 4-7 neurons, we set an initial value of 5 clusters per site. This is a good compromise between nonsupervised clustering and accurate isolation of sources. These results are further compared against "ground truth data" from SIMONE. To illustrate some of the key aspects that make SIMONE a useful tool for validating the algorithms performance, we show: 1) the output of the detection scheme versus the true spiking activity and 2) the output of K-means clustering versus the projection of simulated spikes onto the same PCA space used for detected data.

III. RESULTS

The tests relying on unbiased metrics show a great degree of similarity between real and simulated signals, as shown in Table III. It is remarkable how statistics match since noise has not been directly measured (see Section II-D-2). Electrodes 3 and 4 ( and ) are less similar to the original signal because most of the neighboring neurons are close to the limit of the range of sensitivity. Furthermore, SNR and in those probes are consistent with the decrease of spiking amplitude.

Quantiles are points taken at regular intervals from the cumulative distribution function of a random variable. The quantile-quantile plot (Q-Q plot) is a data analysis technique for comparing the distributions of two sets of data. A 45 reference line is also plotted. If the two sets come from a population with the same distribution, the points should fall approximately along this reference line. This test [see Fig. 6(c)] shows that both datasets are very similar.

Autocorrelation is an important measurement because it gives us information about the intrinsic dynamics of the signal. Both simulated and real data present approximately the same autocorrelation [see Fig. 6(a) and (b)], meaning that their power spectral densities are also alike (the power spectral density of a signal is the Fourier transform of its autocorrelation function). Note that only the autocorrelation is used to estimate the residuals in the predicted regression of the raw signal. Finding the same autocorrelation means that the dynamics of the signals are essentially the same. This is a very encouraging outcome since no information has been provided regarding the topology or spiking regime of the network.

Even though the signals do not have the same amplitude statistics (which can be easily corrected by adjusting the tissue conductivity for instance), the SNR and frequency spectrum are alike. Moreover, we demonstrate that extracellular noise can indeed be modeled by a normal distribution. For that we propose a novel method of estimating noncorrelated noise hidden in extracellular raw signals. We show that even if spikes represent a small portion (30% in our case) of raw extracellular recordings 11 , they have a great impact on the overall dynamics of the signal. Both signals (real and simulated) are presented in Fig. 7. The difference between real and simulated data is mainly due to the stochastic positioning of neurons: we may safely assume that experimental data presents the activity of neurons placed at midrange from electrodes, rather than close-range (as some of those from simulated data). The reader may easily notice that site records some "large" action potentials (these come from neuron 11).

In the following, we present a qualitative description of the SpikeOMatic-based analysis to illustrate how SIMONE can indeed be used as input to spike processing algorithms. Without losing generality, we show results only from site due to space constraints. Following the description of Section II-E, spikes are detected. Detected events during a reduced period of time are marked in Fig. 8(a) along with the true spikes 12 . Note that most of the spikes are accurately detected. The fact that precise spiking activity is known can give the user useful information concerning the performance of the detection process.

Concerning spike sorting, we analyze the clustering obtained using K-means. As stated before, the clustering is given an initial grouping value of 5. The resulting classification may be observed in Fig. 8(c). We use the PCA matrix previously obtained to project the "actual" spikes [Fig. 8(d)]. In this way, it is easy to have a preliminary idea about the classification accuracy. For instance, from this example we observe that there exists a strong correlation between clusters A, B, and E from Fig. 8(c) and clusters linked to neurons 11, 5, and 2 from Fig. 8(d) respectively. Cluster C from the first plot actually covers the activity of three neurons having similar spike templates [START_REF] Smith | A tool for synthesizing spike trains with realistic interference[END_REF][START_REF] Hines | NEURON: A tool for neuroscientists[END_REF][START_REF] Rutishauser | Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo[END_REF]. On the other hand, cluster D is an artifact resulting from superposed spikes from different neurons. Events from plot 8(c) not appearing in Fig. 8(d) correspond to falsely detected spikes. Inversely, those events appearing only on the second plot indicate missed spikes.

IV. CONCLUSION

In this paper, we present a novel tool providing a simple way to simulate realistic MEA-based acquisitions from a fully defined neuronal network. A fundamental requisite is that the number of active cells, their position and the exact firing activity of each of them must be known at all times. Our simulation paradigm is based on two main axes: a coupled model describing the network dynamics and the extracellular transfer model, and a fully statistical parameterization. The dual-layer model is based on commonly used approximations (such as the integrate-and-fire dynamics) and experimental observations. The randomness of the parameters allows for a great number of simulations to be made starting from the same defining functions. This is especially useful when validating the performance of an algorithm under closely related scenari. We explain how the dual-layer model makes this simulation engine suitable for testing the robustness and accuracy of spike processing algorithms (such as the detection, alignment and sorting of spikes). The majority of data that should eventually be treated by a data processing system are simulated: a spike signature for each neuron with time-dependent modulation, the correlated activity between connected neurons, spike burst regimes and realistic noise levels. The system also accounts for acquisition noise and time-dependent degradation of electrode gain in a simple way.

We have validated our tool by accurately reproducing extracellular recordings from the Periplaneta americana antennal lobe acquired in vivo in our facilities. First, we apply several unbiased metrics to compare simulated and real signals. These metrics were specifically chosen to be independent of any a priori information (such as the exact spiking dynamics). Then we show some interesting results that highlight how a spike processing algorithm can be analyzed and tested when the "ground truth data" is perfectly known. We also propose a method to calculate real extracellular noise from the autocorrelation of raw signals.

Naturally, SIMONE has some limitations, such as the inability to quickly simulate large network synchronization, to mimick precise tissue response (e.g., to odors by the antennal lobe), or to reproduce filtering properties of the tissue. But overall, we show that it provides some very interesting characteristics that make it suitable for testing spike processing algorithms and describing the activity in small networks.

Fig. 1 .

 1 Fig. 1. Coupled stages of the dynamic model. (a) Current diagram for the leaky integrate-and-fire model. Conductance models are not shown. (b) Extracellular model in the form of a pseudo monopole of current.

Fig. 2 .

 2 Fig.2. Experimental setup[START_REF] Chaffiol | Quantitative analysis and models of projection neurons responses to odours in the first olfactory relay of an insect, the cockroach Periplaneta americana[END_REF]. Odors are presented through an olfactory delivery system. A main flow of charcoal-filtered and humidified air is continuously delivered through a glass tube (A) in front of the cockroach antennae. Pasteur pipette (B), containing a filter paper with an olfactory stimulus (black arrow), is inserted into the main flow glass tube (A). Stimulus from the Pasteur pipette is controlled by an electrovan.

  Fig. 3. Synaptic kernel . Postsynaptic transient due to a each presynaptic spike vanishes after 25 ms.

Fig. 4 .

 4 Fig. 4. Typical -induced spontaneous spiking for neurons 1-5 during a 2-s time window. Dirac-pulse of current produced each spike event.corresponds to the typical spontaneous spiking frequency for each neuron. Spike waveforms are not shown for simplicity.

Fig. 5 .

 5 Fig. 5. Spike templates generation and typical amplitudes. (a) Experimental extracellular recorded spikes (Fig. 2). (b) Intracellular spike templates given to SI-MONE as an input parameter. These are estimated from waveforms in Fig. 5(a). (c) Simulated extracellular spikes generated by each neuron and "recorded" at a distance of 40 m. (d) Resulting distance-dependent attenuation of the maximum spike amplitude provided by SIMONE for S/m and .

Fig. 6 .

 6 Fig. 6. Comparison tests between real data and simulated signals on electrode . Real data is acquired according to the protocol described in Section II-C. Electrode 1 is arbitrarily chosen (all electrodes showed similar results). (a) Autocorrelograms. (b) Power spectral density. (c) Q-Q plot between simulated data and real data. When datasets have the same distribution the curve is a 45 dark line. (d) Distribution for residuals (bar plot) from the autoregressive model applied on real data. Straight line is the Gaussian distribution ( V; V).

Fig. 7 .

 7 Fig. 7. Compared signals and simulation scenario. (a) Real spontaneous response from the American cockroach antennal lobe. (b) 7 s-long simulated signal. Activity simulated for electrode . Dotted lines in Fig. 7(a) and (b) shows the mean and the threshold used to estimate the SNR. (c) Automatically generated neural network. Neurons are placed in a 2-D plane. MEA is placed at 20 m from the cells. Synapses are not shown. (d) Schematic side view of simulated scenario.

Fig. 8 .

 8 Fig. 8. Spike sorting outcome on simulated data. (a) A 1-s sample from the 10-s simulation. Detected spikes are marked with red dots. Actual events time location is shown below the signal. (b) Extracted spikes (over a 2.5-ms window). They are aligned with their absolute local maxima at 1 ms. (c) PCA projection over the two principal directions. Colors indicate the result of K-means clustering. (d) Actual events projected over the same PCA space as (c). Each neuron is associated with a color and events are colored accordingly.

Ricardo

  Escolá received the degree in electronic engineering from the Universidad Nacional del Sur (Bahía Blanca), Argentina, in 2002. He is currently working toward the Ph.D. degree from the Institut National des Sciences Appliquées de Lyon, France. Since 2003, he has been developing his research activities at CEA/Leti in Grenoble (France) taking part in several multidisplinary projects in neuroscience. His research interests include digital electronics, neural signal processing, and neural interfacing. Christophe Pouzat received the M.Sc. degree in physics and the Ph.D. degree in biology for an electrophysiological investigation of the synaptic connection between interneurons of the molecular layer and Purkinje cells in the cerebellar cortex from Orsay University, in 1992 and 1998, respectively. He then spent two years as a postdoctoral fellow in the laboratory of Dr. G. Laurent at the California Institute of Technology, Pasadena. He has been a research fellow in the Cerebral Physiology Laboratory (CNRS and University Paris-Descartes) since 2001 where he works on MEA recordings from the insect olfactory system and develops neurophysiological data analysis software. Antoine Chaffiol received the M.S. degree in neuroscience, behavior and cognition from Toulouse University, France, in 2001. He is currently working toward the Ph.D. degree in neuroscience in the Laboratoire de Physiologie Cerebrale, CNRS, Paris. His team is studying olfactory representations from neurons in the insect antennal lobe using multielectrode array recordings. Blaise Yvert received the M.Eng. degree in applied physics from the Ecole Centrale de Lyon and Cornell University in 1993, and the M.Sc. and Ph.D. degrees in biomedical engineering from the Institut National des Sciences Appliquées (INSA) de Lyon, in 1994 and 1996, respectively. He is now a Research Fellow at the Center for Integrative and Cognitive Neuroscience in Talence (CNRS and University Bordeaux I & II). His interest focuses on activity-dependent mechanisms underlying the development of neural networks. Isabelle E. Magnin heads the CREATIS-LRMN Laboratory devoted to medical imaging, in Lyon, France. She belongs to numerous french, european and international scientific boards. She is an Honorary Professor of the Harbin University in China. Her research interest is medical image processing including inverse problems, tomography, motion estimation, analysis, and modeling of 2D/3D deformable objects from multimodal image sequences. Since 2000, she promoted biomedical image processing applications in European Healthgrid projects. She is the author of about 230 publications. Dr. Magnin received the gold medal for innovation from the European Society of Non Destructive Testing for original work in ultrasonic image processing, in 1994. Régis Guillemaud graduated from the Institut National Polytechnique de Grenoble, France, in 1985. He was first involved in vision systems for geophysics and for autonomous robots. Since 1991, he has been involved in the field of medical imaging. He developed 3-D reconstruction algorithms (for nuclear medicine and CT) and correction algorithms for MR images. From 1997 to 2002, he was Project Manager for the development of a digital X-ray detector for medical applications. Since 2003, he has been in charge of microsystems-related projects for epilepsy diagnosis and neuroelectronic interfacing.

TABLE II PARAMETERS

 II USED IN THE SPIKE DETECTION/SORTING PROCESS

	TABLE III
	SIGNALS SIMILARITY QUANTIFICATION

bits per second.

This software is not freely available for download at the moment. Prospective users are invited to contact the authors for further information.

See Nomenclature for a glossary of mathematical symbols and notation.

Action potentials fired by the same neuron have little dispersion and may be considered identical[START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF].

Current flow lines are directed radially towards the cell.

NeuroNexus, US: a2x2-tet-5mm150-150-312.

Syntech, NL.

This has been empirically validated during the last couple of years by the Laboratoire de Physiologie Cérébrale (Chaffiol-Pouzat, unpublished data).

Spikes are detected and classified using SpikeOMatic[START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF],[START_REF] Pouzat | Improved spikesorting by modeling firing statistics and burst-dependent spike amplitude attenuation: A markov chain monte carlo approach[END_REF]. Events assigned to the same neuron are averaged to build a single template.

Authors generally assume "worst case classification complexity scenario" of 50 spikes/s per electrode[START_REF] Zumsteg | Power feasibility of implantable digital spike sorting circuits for neural prosthetic systems[END_REF] for human cortical neural activity. This represents roughly 15% of the whole signal.

We only show the first second of simulation but a 10-s simulation is used in the overall process.

ACKNOWLEDGMENT

The authors would like to thank the whole staff involved in the NEUROCOM project for their generous support (A. Bourgerette in particular for the fruitful debates concerning electronic noise); to J. Torbet and J.-F. Bêche for polishing the use of English in the manuscript. The authors would also like to sincerely thank the reviewers for their thorough work providing specific comments, suggestions, and corrections that significantly improved the quality of the paper.