NeuroMap: A Spline-Based Interactive Open-Source Software for Spatiotemporal Mapping of 2D and 3D MEA Data

A major characteristic of neural networks is the complexity of their organization at various spatial scales, from microscopic local circuits to macroscopic brain-scale areas. Understanding how neural information is processed thus entails the ability to study them at multiple scales simultaneously. This is made possible using microelectrodes array (MEA) technology. Indeed, high-density MEAs provide large-scale coverage (several square millimeters) of whole neural structures combined with microscopic resolution (about 50 μm) of unit activity. Yet, current options for spatiotemporal representation of MEA-collected data remain limited. Here we present NeuroMap, a new interactive Matlab-based software for spatiotemporal mapping of MEA data. NeuroMap uses thin plate spline interpolation, which provides several assets with respect to conventional mapping methods used currently. First, any MEA design can be considered, including 2D or 3D, regular or irregular, arrangements of electrodes. Second, spline interpolation allows the estimation of activity across the tissue with local extrema not necessarily at recording sites. Finally, this interpolation approach provides a straightforward analytical estimation of the spatial Laplacian for better current sources localization. In this software, coregistration of 2D MEA data on the anatomy of the neural tissue is made possible by fine matching of anatomical data with electrode positions using rigid-deformation-based correction of anatomical pictures. Overall, NeuroMap provides substantial material for detailed spatiotemporal analysis of MEA data. The package is distributed under GNU General Public License and available at http://sites.google.com/site/neuromapsoftware.

pitch small enough for information to become redundant across neighboring electrodes, the coverage of large neural structures with a limited number of electrodes would benefit from methods that estimate data continuously across the array. For instance, such methods would be of particular interest for mapping local field potentials (LFPs), which have a spatial frequency lower than action potentials and can thus be appropriately interpolated between recording sites.
The issue of spatial interpolation of bioelectric activity has been extensively addressed these past decades in the context of human electrophysiology using electro-and magneto-encephalography (EEG and MEG). Among the different approaches considered for interpolation, spline interpolation was found to be one of the most attractive compared to other interpolation methods (Perrin et al., 1987;Yvert et al., 2001). Indeed, spline interpolation has several properties of major interest: (1) it does not require equally spaced spatial sampling and can thus be used for any arrangement of the recording sites; (2) it does not constrain interpolated extrema to the range of sampled values, unlike bilinear and nearest-neighbors interpolations; (3) it is computationally simple, requiring no more than the resolution of a linear system of equations (Harder and Desmarais, 1972); (4) it provides an analytical expression of the interpolant, which is differentiable and convergent (Duchon, 1976). Differentiability of the interpolant is interesting for the computation of the spatial Laplacian, allowing to estimate the current source

IntroductIon
Microelectrode arrays (MEAs) are a powerful tool to record large populations of neurons and monitor local neural networks as a whole. They are increasingly employed in a wide variety of studies, in vivo (Nicolelis et al., 1997;Csicsvari et al., 2003) as well as in vitro (Shahaf and Marom, 2001;Tscherter et al., 2001;Beggs and Plenz, 2003;Wagenaar et al., 2005;Steidl et al., 2006;Rolston et al., 2007). Progress in the field of microtechnologies now allows the development of high-density arrays, opening the route toward exhaustive spatial sampling and coverage of the explored neural structures (Berdondini et al., 2005;Charvet et al., 2010). In particular, systems handling high number of channels now offer the possibility to sample the neural tissue in 3D with noticeable numbers of recording sites along all dimensions in space (Perlin and Wise, 2008;Rousseau et al., 2009). Yet, few tools are available to take advantage of highdensity 2D or 3D MEA data to provide synthetic spatiotemporal representations of activity across all the channels, or to localize the underlying active areas on the tissue anatomy.
Synthetic representation of raw potential data or processed neural signals (e.g., spike firing rates) across all the channels of an array can be achieved by building images of activity at a given time point. Classically, activity maps are built using as many pixels as electrodes. However, this approach provides scattered information at recording sites only, and does not provide any estimate of activity between electrode sites. While these methods are satisfactory for electrode A major characteristic of neural networks is the complexity of their organization at various spatial scales, from microscopic local circuits to macroscopic brain-scale areas. Understanding how neural information is processed thus entails the ability to study them at multiple scales simultaneously. This is made possible using microelectrodes array (MEA) technology. Indeed, high-density MEAs provide large-scale coverage (several square millimeters) of whole neural structures combined with microscopic resolution (about 50 μm) of unit activity. Yet, current options for spatiotemporal representation of MEA-collected data remain limited. Here we present NeuroMap, a new interactive Matlab-based software for spatiotemporal mapping of MEA data. NeuroMap uses thin plate spline interpolation, which provides several assets with respect to conventional mapping methods used currently. First, any MEA design can be considered, including 2D or 3D, regular or irregular, arrangements of electrodes. Second, spline interpolation allows the estimation of activity across the tissue with local extrema not necessarily at recording sites. Finally, this interpolation approach provides a straightforward analytical estimation of the spatial Laplacian for better current sources localization. In this software, coregistration of 2D MEA data on the anatomy of the neural tissue is made possible by fine matching of anatomical data with electrode positions using rigid-deformation-based correction of anatomical pictures. Overall, NeuroMap provides substantial material for detailed spatiotemporal analysis of MEA data. The package is distributed under GNU General Public License and available at http://sites. google.com/site/neuromapsoftware. densities (CSDs) of the extracellular potential field (Mitzdorf, 1985;Nunez and Pilgreen, 1991;Law et al., 1993;Nunez et al., 1994;Babiloni et al., 1996).
Because such approach has, to our knowledge, not been adapted for mapping MEA data, here we present NeuroMap, an interactive Matlab-based software for mapping MEA data based on spline interpolation. NeuroMap presents itself in a graphical, user-friendly, interface allowing complete parameterization of algorithms, and visual renderings. It can handle data from any MEA configuration, including irregularly spaced electrodes arrays with electrode positions distributed in 2D or 3D. In addition to potential mapping, Neuromap also offers the possibility to map spatial Laplacian of data to provide CSD estimates. Another key feature in the case of 2D data is the possibility to coregister activity data onto anatomical images. NeuroMap has been used on real LFP data recorded from embryonic mouse brainstem and spinal cord preparations. NeuroMap is made freely available to the community under GNU General Public License (GPL) at http://sites.google.com/site/neuromapsoftware.

MaterIals and Methods 2d surface splIne InterpolatIon
Let n be the number of recording sites and (x i , y i ) 1≤i≤n their coordinates. The estimate of the potential at any given location (x, y) is given by the following mth degree surface spline expression (Duchon, 1976;Perrin et al., 1987): where k m is the function defined by: If (V i ) 1≤i≤n is the vector of the potentials at the electrodes, Q = (q 00 , q 01 , q 11 ,…, q m−1 , m−1 ) T and P = (p i ) 1≤i≤n , then P and Q are obtained through the resolution of the following system of linear equations: where ...

3d voluMe splIne InterpolatIon
Spline interpolation can be generalized in higher dimensions by means of minor changes (Duchon, 1976). In three dimensions, the interpolated spline function becomes (Babiloni et al., 2000;Yvert et al., 2001): where h m is the function defined by: f m depends on the n coefficients p i (vector P) and the m(m + 1) (m + 2)/6 coefficients q dkg (vector Q). P and Q are obtained through the resolution of the following system:

resolutIon and regularIzatIon
In both 2D and 3D cases, spline interpolation requires resolution of a system of linear equations (Eqs 2 or 4) of the type AX = B. A numerical approximate solution  X of this system is: where A + is the pseudo-inverse of A, classically obtained using the singular value decomposition of A: with Σ + a diagonal matrix containing k non-zero values being the inverse of the k highest singular values of A. In most cases, the value of k is the number of non-zero singular values of A. However, when A is highly ill-conditioned, a smaller value of k can be considered to obtain a regularized solution to the linear system (Uutela et al., 1999), an option available in NeuroMap.

spatIal laplacIan
As spline functions of mth degree have well-defined and continuous m − 1 derivatives, spatial Laplacians are readily available for spline interpolants of third degree or higher. In both surface (2D) and volume (3D) interpolations, there are two ways of computing the spatial Laplacian. One may compute the values of the Laplacian on the grid points, using the analytical expression of the spline function Laplacian (derived from Eqs 1 or 3): 1≤i≤N are the coordinates of the N grid points.
Another way is to compute a discrete and numerically approximated Laplacian of interpolated data: respectively. The potential values at the electrodes were computed from current sources as indicated above, such as ∆ cs = 2.5∆ e and Z cs = −∆ e /2(see Figure 1, top row). Likewise, 3D fictive sets of data were created from virtual 3D MEAs of 125 (5 × 5 × 5) or 729 (9 × 9 × 9) equally spaced electrodes, and with inter-electrode distances ∆ e and ∆ ′ e respectively. The electrode located at the center of these arrays was placed at the origin of the (x, y, z) coordinate system. The potential values were computed from the same sources (x 1 = −1.25∆ e , y 1 = 0, z 1 = −∆ e /2 and x 2 = 1.25∆ e , y 2 = 0, z 2 = −∆ e /2; see Figure 1, bottom row).
For our simulations we used ∆ e = 200 μm, but defining ∆ cs as a function of ∆ e makes our results scale-independent.

separatIon power
The possibility of improving source localization using the spatial Laplacian was investigated from the simulated data described above. The two parallel current sources create two potential fields, each one reaching its maximum amplitude in a distinct hemispace. Viewed from the recording sites, and even after careful interpolation, these two potential fields might be more or less easily differentiated, depending on the separation between the two sources (∆ cs ) in relation to the electrode pitch (∆ e ). In order to estimate the quality of this discrimination when interpolating the potential field versus interpolating the spatial Laplacian of the potential field, we define an index of separation power (SP) as follows: where V max is the maximum amplitude reached in each hemispace and V mid the value at the point located between the two local maxima (see Figure 5C). The higher SP, the better the separation of the sources inferred from the interpolated map.
OF1 mice (Charles River Laboratories, L'Arbresle, France). As previously described (Joucla and Yvert, 2009), the neural tube was opened dorsally and meninges were removed prior to positioning on MEAs. Experimental protocols conformed to recommendations of the European Community Council and NIH Guidelines for care and use of laboratory animals. Microelectrode array recordings were performed using two different systems. First, LFPs were recorded using a rectangular 4 × 15-electrode array chip from Ayanda BioSystems Inc. (Lausanne, Switzerland), amplified using MultiChannelSystems MEA1060 amplifier (MCS, Reutlingen, Germany), and digitized at 10 kHz using CED Power1401 and Spike2 software version 6 from Cambridge Electronic Design (CED, Cambridge, UK). Second, a non-rectangular 256-electrode array fitting the dimension of the embryonic hindbrain and made by ESIEE-Paris was connected to the BioMEA system (Charvet et al., 2010) now commercialized by Bio-Logic SAS (Grenoble, France) and signals were digitized at 13 kHz. In all cases, data were low-pass filtered below 50 Hz offline.

sIMulated 2d and 3d data
Simulated data consisted in the potential field created by two current dipoles of identical moment and orientation located at positions (x 1 = −∆ cs /2, y 1 = 0, z 1 = Z cs ) and (x 2 = ∆ cs /2, y 2 = 0, z 2 = Z cs ). Their moments p 1 and p 2 were of unit strength oriented along the positive z-axis (p 1 = p 2 = u z ). The potential values at locations (x, y, z) were computed as in an infinite homogeneous volume conductor using the following relation (Malmivuo and Plonsey, 1995): Two-dimensional fictive sets of data were created from virtual 5 × 5 and 9 × 9 planar MEAs located in the y = 0 plane, centered on (x = 0, z = 0), and with inter-electrode distances ∆ e and ∆ = ∆ ′ e e /2 Figure 1 | 2D and 3D maps of fictive MeA data. Left: current sources positions and orientations as described in Methods. Middle: simulated data computed on grids of 1681 and 68921 points for 2D and 3D configurations respectively. Right: data from fictive 5 × 5, 9 × 9, 5 × 5 × 5 and 9 × 9 × 9 MEAs interpolated on 1681 (2D) or 68921 (3D) points grids and mapped by NeuroMap. The potential field geometry can be estimated by interpolation. Note that the interpolated map is best for highest MEA resolution. Due to the strong dynamic aspect of bioelectrical signals, it is essential to explore data also along the time dimension. A classical representation of both spatial and temporal features of electrical activity is time curves. NeuroMap provides such a representation, where time curves of all recording sites are positioned according to the MEA geometry ( Figure 3A). It has the advantage of allowing immediate spotting of defective electrodes, which can then be disabled. Data from disabled electrodes are not used in mapping interpolation.

IMage warpIng
For coregistration of activity onto anatomical pictures, fine image warping is sometimes required for accurate matching of picture with recording sites locations. For this purpose, an image warping tool was implemented in NeuroMap, based on a moving least squares algorithm employing rigid transformations for a realistic deformation (Schaefer et al., 2006). This warping strategy consists in defining a set of control points p and their deformed positions q. The algorithm then constructs a deformation function f defined as This warping tool works in an interactive fashion, allowing the user to control the deformation through a set of handles that can be freely defined and positioned by mouse (see Figure 6).

the neuroMap package
NeuroMap is licensed under GNU GPL and can be freely downloaded from a dedicated website (http://sites.google.com/site/ neuromapsoftware). NeuroMap comes as a graphical user interface developed under Matlab and compatible with versions 7.1 and higher. NeuroMap has its own data format which specifies the names and 3D coordinates of the recording sites, followed by data recorded along time at each recording site. We provide NeuroMap with routines to write NeuroMap data files from Spike2 data files or Matlab variables. Besides, the distribution package includes data sets used in this paper and a comprehensive documentation.

results
Examples of 2D and 3D maps of fictive MEA data generated by two current dipoles (see Materials and Methods) are provided in Figure 1 together with the exact calculated potential field. It appears that even with a number of electrodes as small as 25, a good approximation of the potential field pattern can be obtained using spline interpolation.
The software includes a graphical user interface (Figure 2A) allowing complete configuration of data mapping: selection of single or multiple times of mapping; selection of mapping mode (raw data, analytical, or numerical Laplacian); parameterization of the spline interpolation; loading and matching of anatomical picture with electrodes; configuration of visual rendering, including transparency effects, colormap aspect, and miscellaneous plottingrelated parameters.
Examples of 2D mapping using experimental data from a 60-electrode rectangular MEA and from a 256-electrode nonrectangular MEA are provided in Figures 2B,C, respectively. A spatiotemporal representation of the data consists in building a series of maps at regular time intervals (Figure 3B). Such a representation is more convenient to visually detect various spatiotemporal features (initiation site, directions of propagations, etc.) than time curves.

voluMe representatIon of 3d Mea data
Three-dimensional sampling of neural tissues is promising to more accurately localize activity in neural structures, as for instance to discriminate between cortical layers involved in different types of behaviors (Krupa et al., 2004). Upcoming MEA systems promise to make it a reality (Perlin and Wise, 2008;Rousseau et al., 2009). Anticipating the emergence of 3D-recording technologies, NeuroMap was designed to offer a user-friendly interface for exploration of volume data (Figure 4). All configuration possibilities offered in 2D are available in 3D, including temporal and spatiotemporal representation of data in 3D views or any x, y, or z planes in the volume of the array (not shown here).

spatIal laplacIan
Spline Laplacian was shown in various human electrophysiological studies to increase significantly the spatial resolution of EEG potential distributions by providing an estimate of the cortical surface potential (Nunez and Pilgreen, 1991;Law et al., 1993;Nunez et al., 1994;Babiloni et al., 2000). In addition, providing some assumptions on the medium conductivity, one may show that the spatial Laplacian of the extracellular potential is directly linked to the current source density (Mitzdorf, 1985). Here, we implemented 2D and 3D spatial Laplacian of MEA data, which is illustrated in Figure 5 in the case of 3D simulated data (cubic 5 × 5 × 5 MEA). Raw potential and spatial Laplacian data were compared for their ability to separate the two current sources used to simulate the potential field.
For distances between sources comparable to electrode pitch, the local potential extrema generated by the parallel current sources appear merged (Figure 5A, top-left), but they could be distinguished using spatial Laplacian (Figure 5A, bottom-left). When distance between sources increased, the map of raw potential could hardly display two distinct local maxima in a unique field (Figure 5A, top-right), while the Laplacian performed complete, clear-cut separation (Figure 5A, bottom-right). For large distances between sources, the spline Laplacian tended to generate phantom activity between the two fields ( Figure 5B, top) due to rebounds of the spline, a phenomenon that could be corrected by regularization (Figure 5B, bottom; see Resolution and Regularization in Materials and Methods). data. NeuroMap offers a convenient way to do so by coregistration of maps of activity onto images, which are typically pictures of the recorded preparation ( Figure 6A1). To perform this coregistration, the software requires only information about the exact localization of the electrodes on the picture. This is determined from the positions of two reference electrodes, selected and indicated on the picture by the user (Figures 6A2,A3).
When large structures are considered, the full image is usually obtained by stitching several smaller fields of views, a process that may introduce distortions in the final anatomical image reconstruction. In order to compensate for these distortions and achieve accurate registration of data mapping on the anatomy of the tissue, image warping can be performed (Figures 6A4,A5).
The performance of source separation was further investigated systematically for intersource distance ∆ cs ranging from 0 to 4 ∆ e by 0.05 × ∆ e steps. Figure 5C illustrates how the SP index was calculated (see also Materials and Methods). Figure 5D displays SP values in the whole ∆ cs range for potential and Laplacian interpolations. Spatial Laplacian always improved sources separation when compared to raw potential. While regularized Laplacian displayed a lower SP than original Laplacian, it remained more effective than raw potential for source separation.

coregIstratIon on anatoMIcal IMages
A major feature of the MEA technology is the possibility to cover a large network and to study interactions between various areas. For this purpose, it can be useful to merge signal data with anatomical To each view is associated a set of controls (a) for printing, temporal mapping, and video making. Extra controls (b) come with the volume view to modify the volume display mode: orthogonal slices (as shown here) or parallel slices; in the latter, direction and number of slices can be freely specified.

dIscussIon
We presented NeuroMap, an interactive software for spatiotemporal mapping of 2D and 3D MEA data. NeuroMap can be fed with any quantity at each electrode, which can be either raw bioelectric data or processed data derived from recordings (e.g., spike firing rates, LFPs, etc.). To build maps, NeuroMap estimates data between recording sites using thin Matching electrode locations onto anatomical picture is possible for any geometry of MEA. An example of 256-channel activity coregistered onto the anatomy of an embryonic mouse hindbrain is shown in Figure 6B. The visual rendering of coregistration can be further improved by applying transparency to activity data (not shown here), or by thresholding mapped values so as to emphasize regions of highest activity, as in Figure 6B. Top: map illustrating the phantom activity that can appear in the spatial Laplacian (red arrowheads). Bottom: phantom activity vanishes after regularization. (C) Illustration of how the SP index is calculated. In the simulated data, parallel current sources generate maximum potential field in the y = 0 plane (top). The line on which absolute maximum is reached is used. The SP index is defined using the local extrema on the maximum line (bottom). It has a value of 0 when the potential along the line has only one maximum, and is close to 1 for an optimal separation of the maxima. (D) Evolution of SP according to distance between sources (∆ cs ) for raw potential (blue), spatial Laplacian (red), and regularized Laplacian (orange). Circles indicate maps drawn in (A,B) layouts that may be developed for specific uses or neural structures (e.g., hippocampus, brainstem, etc.). It is as well adapted to those forthcoming 3D arrays that will provide recording sites in all directions of space. Spatial Laplacian is an effective approach to further increase spatial resolution and accurately estimate the locations of the current sources that underlie the observed extracellular potential plate spline interpolation. This method provides a smooth estimation of data and can predict the presence of extrema between the electrodes, unlike nearest-neighbor or bilinear interpolation approaches. The use of this software is not restricted to classical rectangular and planar commercial systems, nor to a regular spacing between electrodes. NeuroMap was designed to deal with any geometry of arrays, so as to handle the whole variety of MEA The picture of a whole embryonic hindbrain-spinal cord preparation recorded with a 4 × 15 MEA was rebuilt from six partial pictures (dashed frames). As partial pictures have not been perfectly aligned during reconstruction, the anatomical picture needs appropriate warping to match the real MEA geometry. (A2) First, user identifies two or more electrodes on the picture (white crosses). (A3) Using the identified electrodes coordinates, the program computes the scale of the picture and the theoretical positions of all electrodes (white crosses). As expected, these do not match the electrodes of the picture. (A4) In order to perform appropriate deformation of the picture, the user specifies the correspondence between theoretical and pictured positions for a few electrodes (white circles). These landmarks determine reference pixel displacements that are extrapolated for the whole image. (A5) Result after application of the deformation algorithm. (B) Mapping of activity on a matched anatomical picture of an embryonic hindbrain (data is identical to that of Figure 3B). In order to improve visibility of the anatomy, values in the noise range (±4 μV) were not displayed.