Lorraine Perronnet

Maria Eugenia Vilarchao

Guillaume Hucher

Daniel E Shulz

Gabriel Peyré

Isabelle Ferezou
email: isabelleferezou@unic.cnrs-gif.fr

María Eugenia Vilarchao

An automated workflow for the anatomo-functional mapping of the barrel cortex

Keywords: barrel cortex, histological sections, blood vessels, registration, robust iterative closest point, gradient domain fusion, voltage sensitive dye imaging

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The rodent primary somatosensory cortex is a very convenient model for studying the cortical processing of sensory information because of its well defined structural and functional layout that is invariant from animal to animal [START_REF] Welker | Quantitative correlation between barrel-field size and the sensory innervation of the whiskerpad: a comparative study in six strains of mice bred for different patterns of mystacial vibrissae[END_REF][START_REF] Meyer | Cellular organization of cortical barrel columns is whisker-specific[END_REF][START_REF] Egger | 3D reconstruction and standardization of the rat vibrissal cortex for precise registration of single neuron morphology[END_REF]. In its layer IV, neurons are gathered into clusters called barrels that respect the same topology as the whiskers on the snout of the animal [START_REF] Woosley | The structural organization of layer IV in the somatosensory region (SI) of mouse cerebral cortex: The description of a cortical field composed of discrete cytoarchitectonic units[END_REF]. Each barrel is dedicated primarily to the processing of the input coming from its corresponding whisker (Figure 1A,B). When studying sensory processing in the barrel cortex either with electrophysiological or imaging methods, it is therefore of great interest to superimpose the recorded activ-ity onto the underlying barrel topography, which can be reconstructed from the post hoc alignment of tangential brain slices stained for cytochrome oxidase. In order to optimize this anatomo-functional mapping, which is usually accomplished manually, we developed an automated workflow for the registration of the histological slices of the barrel cortex and the 2-D barrel map reconstruction.

Here we focus our attention on voltage sensitive dye imaging (VSDI) of the mouse barrel cortex to illustrate the usefulness of the approach. However, the method can be extended to the study of the rat barrel cortex and applied to other techniques such as 2-photon calcium imaging.

The traditional first step to recover the map of the barrel cortex after imaging experiments is: brain fixation by perfusing the animal with a solution of paraformaldehyde, followed by the cutting of tangential slices (∼100 µm thick, with or without previous flattening of the cortex), which are subsequently stained for cytochrome oxidase using classical histological procedures that reveal the barrel arrangement in layer IV (Figure 1C,D [START_REF] Land | Cytochrome oxidase staining in the rat smi barrel cortex[END_REF]). A. Following the registration of tangential histological slices and the reconstruction of the barrel map, one can see that the spatial organization of the layer 4 barrels matches the layout of the vibrissae on the snout of the animal (B). C. Drawing of a coronal section of the left hemisphere of the mouse brain illustrating the position of layer 4 barrels within the primary somatosensory area of the cortex. After in vivo imaging of barrel cortex activity, sections are cut tangentially to reconstruct the layer 4 barrel map (cutting plane indicated by the red line). D. A series of tangential histological slices stained for cytochrome oxidase. On the first slice one can see superficial blood vessels. On the other slices, one can see white circular to elliptic spots that correspond to sections of plunging blood vessels. Depending on the exact axis of the cut, barrels can be spread over several slices.

Next, using digital microphotographs of the slices, it is necessary to:

1. register the slices ; 2. fuse the registered slices to define a reconstructed barrel image.

In this article we provide an automated solution for these two steps which are the most time-consuming tasks of the workflow when using conventional manual methods. After completing these steps, it is then relatively simple to define the barrel map by segmenting the reconstructed barrel cortex image. The superimposition of the map with the functional data can be finally achieved by using the superficial blood vessels as anatomical landmarks (Figure 2). The proposed anatomo-functional mapping tool significantly speeds up the overall process and provides more accurate anatomo-functional mapping.

Registration of histological slices

A typical example of a series of images obtained after the histological process is shown in Figure 1D. Depending on their depth, the histological sections present different properties: in the first section, which corresponds to the surface of the cortex, some large superficial blood vessels are visible together with plunging blood vessels. This section is crucial for the whole process since it contains most of the superficial blood vessels that will be used for the final alignment of the histological data with the VSDI data (Figure 2). The intermediate sections usually contain orthogonal blood vessels (white dots on the slices shown in Figure 1D) that can be used to align subsequent sections. Barrels start appearing on section 3 or 4 and can be visualized on up to 5 sections.

We do not use the deeper sections from subgranular layers as they do not contain any useful information for the reconstruction.

Image registration is a classical problem and its solutions find many applications in medical image analysis. Depending on the imaging modality and the specific prior knowledge of the object to register, a wide variety of methods have been considered in the literature (for an overview see [START_REF] Glocker | Deformable medical image registration: Setting the state of the art with discrete methods[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]). Only a few previous studies explicitly deal with the problem of registering rodent brain histological sections, usually in order to reconstruct a whole brain in 3D [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF][START_REF] Ju | 3d volume reconstruction of a mouse brain from histological sections using warp filtering[END_REF].

In order to exploit the microscopic-scale information of the histological data, these applications require the precise registration of a large number of histological slices from the whole brain with a method that accounts for global but also local nonlinear deformations due to tissue shrinkage and tearing after histological preparation. Ourselin et al. presented in [8] a block matching strategy to compute local similarities and then estimate the rigid transformation that matches the maximum of similar regions in a robust way. Alternatively Ju et al. [START_REF] Ju | 3d volume reconstruction of a mouse brain from histological sections using warp filtering[END_REF] used a method based on pairwise elastic imaging warps, with the specificity to compute the deformation of each section by considering not only two neighbors for each section, but an extended neighborhood including a group of images.

The specific problem of registering histological sections of A. Barrel map reconstructed from the histological slices shown in figure 1D, using our anatomo-functional mapping method. B. Superimposition of the registered histological slices, allowing the delineation of the superficial blood vessels (in blue). C. The superficial blood vessels also appear in vivo on the fluorescent images taken during the VSDI session. The VSDI data can therefore be aligned with the underlying structural map of the barrel cortex using the blood vessels as anatomical landmarks. On the right, the cortical activity is shown imaged at 10ms following a single C2 whisker deflection with the voltage sensitive dye RH1691 under urethane anesthesia. The barrels outlined from the reconstruction in A are superimposed on the image as white lines.

the rodent barrel cortex has been rarely addressed in the literature. Egger et al. [START_REF] Egger | 3D reconstruction and standardization of the rat vibrissal cortex for precise registration of single neuron morphology[END_REF] proposed a tool for the 3-D reconstruction and standardization of the rat barrel cortex for the precise registration of single neuron morphology. Seventy micrometer thick tangential sections of rat barrel cortex are aligned pairwise by finding the rotation and translation that best superimpose blood vessels of two adjacent sections either manually or automatically, using a tool originally developed for the reconstruction of neuronal processes [START_REF] Dercksen | Automatic alignment of stacks of filament data[END_REF]. Finally the barrels are segmented by using a semi-automated method.

Here we developed a tool to compute the rigid transformations to align sets of detected blood vessels (Figure 3A), and we decided to focus on a 2-D image fusion using an automated method. The usual approach to perform point cloud registration is the iterative closest point (ICP) algorithm introduced by Besl [START_REF] Besl | A method for registration of 3-D shapes[END_REF]. While the initial formulation is not robust to outliers, several approaches explicitly deal with this issue [START_REF] Chetverikov | Robust euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF][START_REF] Nishino | Robust simultaneous registration of multiple range images[END_REF][START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Kaneko | Robust matching of 3D contours using iterative closest point algorithm improved by m-estimation[END_REF][START_REF] Ma | Spotlights: A robust method for surfacebased registration in orthopedic surgery[END_REF]. These approaches are related to reweighting least squares methods and we propose in Appendix B a unifying presentation and convergence analysis of a robust ICP method.

Fusion of histological slices

In our framework, the 2-D reconstruction of the barrel maps amounts to perform a fusion of the registered sections (Figure 3B). The goal is to reconstruct the edges of the barrels that are spread over the slices. Image fusion is classically addressed in the field of image processing and computer vision, to perform for instance image editing and stitching. To reconstruct sharp edges, it is necessary to use non-linear methods, the most popular one being based on gradient-domain blending [START_REF] Pérez | Poisson image editing[END_REF][START_REF] Raskar | Image fusion for context enhancement and video surrealism[END_REF]. The goal is to generate a novel image by locally keeping the content of the image having the highest local frequency. In this article, we focus on the use of gradient-domain methods, which have the advantage of being simpler and can be easily tuned for the purpose of barrel map fusion.

Contributions

Our main contribution is a comprehensive pipeline for the 2-D reconstruction of the barrel cortex from tangential histological sections. This pipeline is composed of 2 successive modules that respectively perform: histological section registration and barrel image reconstruction using image fusion. As a side contribution, we propose to recast the registration problem using a robust ICP optimization method, and we show that a Majorize-Minimize framework can be applied to provably ensure the convergence of the method. These contributions originate from specific needs raised by studies of sensory integration in the barrel cortex, however the histological section registration tool proposed here might be helpful to reconstruct any anatomical tissue in which blood vessels penetrate predominantly orthogonally to the cutting plane of the histological slices.

Material and methods

Overview of the proposed framework

This Section details each step of the proposed method. The corresponding Matlab code to reproduce the figures of this article is available online1 . This code is packaged as a graphical user interface (GUI) that is helpful to guide the user through the various processing steps, from image loading to the final reconstructed barrel map.

Table 1 reviews the notation introduced in the paper. Apply the transforms to obtain the registered images { S 1 , . . . , S Q } (Equation ( 5)).

For each m, inpaint the detected vessel and denoise the resulting images to obtain { S 1 , . . . , S Q } (Section 2.3.1). Fuse the relevant maps { S m } m∈M (where M ⊂ {1, . . . , Q} indexes the maps containing apparent barrels) to obtain the fused map S (Section 2.3.2).

Remove the low frequency drift to obtain S (Section 2.3.3).

Step 1: Pairwise Registration of Histological Sections

The input of the algorithm is Q raw images which are digital images with range normalized in [0, 1] (0 being black and 1 white). Figure 1D shows examples of such images.

Pre-processing

Depending on the configuration of the microscope used to acquire the images, the input images contain 2 or 3 regions:

In the center of the image, a gray region corresponding to the histological section of the cortex. This is the foreground region.

A circular white region corresponding to the lens/plate of the microscope and that surrounds the cortex region. It is considered as background.

In some cases, a black region can also surround the two other regions. It is also considered as background. We extract the background using a K-means algorithm with K = 3 regions. The background is then replaced by the value 0 and the resulting images (the so-called "slices" in the following) are denoted {S 1 , . . . , S Q }.

Vessel Detection

Blood vessels which are approximately orthogonal to the cutting plane look like slightly elliptic spots that can be approxi-mated by small Gaussians. To be invariant to local contrast fluctuations, we detect these orthogonal vessels using normalized cross correlations with Gaussian templates of varying standard deviations, assumed to be smaller than σ max .

For each slice S m , for m ∈ {1, . . . , Q}, we compute its associated normalized cross correlation NCC(S m ) against the set of Gaussian templates, as detailed in Appendix A. Given a threshold τ ∈ [0, 1] and a maximum number N of detected vessels, we define the detected vessel centers X m to be the set of pixels x satisfying both NCC(S m )(x) > τ and NCC(S m )(x) is among the N largest values of NCC(S m ).

Slice Registration by Robust ICP

For each m, we now register the slice S m with the slice S m+1 (see Figure 3A). Registration is obtained by computing an optimal transformation T m which maps pixels in slice S m+1 to pixels in image S m . We restrict here the computation to rigid transformations, i.e. of the form T (x) = R(x) + t where R is a planar rotation and t ∈ R 2 is a translation vector.

Variational registration. This optimal deformation is obtained by exploiting the fact that detected orthogonal vessels should have approximately the same position in two consecutive frames. We denote X m+1 = {x i } i∈I and X m = {y j } j∈J the two sets of vessel positions. We cast this problem as the optimization of a non-convex functional measure of the goodness of fit between the transformation T (x i ) of each detected vessel x i in S m and its closest neighbor y j in S m+1 . T m is obtained by computing a local minimizer of min

T i∈I min j∈J ρ(||T (x i ) -y j ||). (1)
Here ρ : R + → R + is a penalty function. The most common choice is a quadratic loss ρ(r) = r 2 , which assumes some sort of Gaussian distribution of the fitting errors. This choice poorly handles outliers in the detected vessels, which are likely to be present in our datasets. We choose here to use the following robust loss function

ρ(r) = log(ε 2 + r 2 ), (2)
which gives less weight to outliers (large values of r) than a quadratic loss. Small values of ε are used to cope with many outliers. Note that setting ε → +∞ recovers the quadratic loss r 2 which assumes no outliers. Note also that other loss functions could be used as well, as long as they satisfy the hypotheses exposed in Appendix B.

ICP iterations. A classical algorithm to minimize (1) is the Iterative Closest Point (ICP), introduced by [START_REF] Besl | A method for registration of 3-D shapes[END_REF] for the quadratic loss ρ(r) = r 2 . This algorithm has been extended by several authors to cope with robust loss (see Section 1.1 for more details). We use a similar approach here, and provide more details in Appendix B.

The ICP algorithm iterates between two steps. In the first step, T is known and assumed to be fixed, and one computes a A. Snapshots extracted from the GUI illustrating the results of the automated detection of blood vessel cross sections on consecutive histological slices (yellow filled circles and red open circles, respectively -S 3 to S 6 previously shown in Figure 1D). The result of the registration obtained by robust ICP applied on the detected blood vessels is illustrated on the right for each pair of slices. Note that the user can delimit manually a region of interest (black dashed lines) in order to restrict blood vessel detection to the barrel field region. B. Images from the same experiment obtained by using our automated barrel map reconstruction tool before (left) or after (right) post-processing steps.

nearest neighbor z i = y j for each vessel x i , where the index j minimizes min

j∈J ρ(||T (x i ) -y j ||). (3)
In the second step, the optimal T is updated by solving min

T i∈I ρ(||T (x i ) -z i ||). (4)
For the quadratic loss ρ(r) = r 2 , this second step is solved in closed form as detailed in Appendix B.2. For a generic loss ρ, there is no such closed form. We detail in Appendix B.1 a Majorize-Minimize (MM) method to compute a local minimizer. To the best of our knowledge, this presentation, and the corresponding convergence analysis, is new.

Initialization. A major difficulty to solve ( 2) is that it is highly non-convex, and the ICP algorithm is likely to converge to a local minimizer T . To improve the quality of the result, it is important to test several initializations to obtain a good registration.

Registered images. Once the registration transforms {T 1 , . . . , T Q-1 } have been computed, they can be cascaded to warp the input slice images to obtain the sequence { S 1 , . . . , S } of sections, all registered with respect to the initial one S 1 = S 1 as follow

S m (x) = S m (T m (x)) where T m = T 1 • . . . • T m-2 • T m-1 . (5)

Step 2: reconstruction of the barrel image

We now have a set { S 1 , . . . , S Q } of registered slices. We fuse them in a single image S which gathers the edge information of the relevant images to reconstruct the barrel map.

Pre-processing

In order to avoid the amplification of artifacts during the gradient fusion process detailed next, we inpaint (i.e. remove) the orthogonal vessel traces and denoise the resulting image. The inpainting method is detailed in Appendix C. The output of the inpainting is then denoised using a median filter on 3 × 3 patches. We use this non-linear filter to reduce salt-and-pepper noise instead of a convolution, as it removes noise while preserving edges. We denote the output of this pre-processing S m .

Slice gradient fusion

We denote m ∈ M the set of relevant slices containing partial barrel information. To reconstruct a sharp image, we fuse together the gradient of the input slices { S m } m∈M by keeping at each pixel the gradient with the largest magnitude. This method is partly inspired by some recent works in computational photography, such as [START_REF] Pérez | Poisson image editing[END_REF][START_REF] Raskar | Image fusion for context enhancement and video surrealism[END_REF]. The details of this method are given in Appendix D. We denote S the output of the fusion process.

Drift removal

The histological sections often present variations in intensity across the barrel cortex that might be due to anatomical reasons. Usually, the anterior lateral barrel subfield (small barrels corresponding to small vibrissae) appears darker than the posterior medial barrel subfield (large barrels corresponding to large vibrissae). This drift is enhanced by the gradient fusion operation that is applied on each pair. As a consequence, the merged image S obtained by the procedure explained above exhibits a strong drift in intensity. We thus filter the merged image with a high-pass filter to remove this low frequency component S = S -S h σ [START_REF] Glocker | Deformable medical image registration: Setting the state of the art with discrete methods[END_REF] where h σ is a low-frequency gaussian filter of standard deviation σ , and is the discrete 2-D convolution. Whisker stimulation. Individual deflections of the right 24 posterior macrovibrissae of the mice were performed using a custom built multi-whisker stimulator based on a matrix of 24 multidirectional piezoelectric benders [START_REF] Jacob | The matrix: A new tool for probing the whisker-to-barrel system with natural stimuli[END_REF]. The whiskers were inserted in 27G stainless steel tubes attached to the benders, leaving 2 mm between the tip of the tube and the whisker base. The 24 whiskers were stimulated individually, in the 4 cardinal directions, at 0.1 Hz within pseudo randomized sequences containing extra blank trials (each stimulation being repeated 10 times). Each whisker deflection consisted of a 100 µm displacement (measured at the tip of the tube), with a 2 ms rising time, a 2 ms plateau and a 2 ms fall (specific filters were used to correct for the mechanical ringing of the stimulators).

Image analysis. Acquisition and data preprocessing were done using in-house software (Elphy, G. Sadoc, UNIC-CNRS), further analyses were made using custom written routines in Ig-orPro (Wavemetrics). Subtraction of the averaged unstimulated blank trials was used to correct for bleaching artifacts. For each whisker, data corresponding to the 4 directions of deflection were averaged.

DiI in vivo staining. DiI stain (Molecular Probes, LifeTechnologies) was deposited on the shanks of silicon electrodes [START_REF] Dicarlo | Marking microelectrode penetrations with fluorescent dyes[END_REF] that were inserted in the barrel cortex perpendicularly with a microcontroller (Luigs & Neumann).

Histological Procedures

Following the experiments and the administration of an overdose of urethane, mice were perfused with saline followed by paraformaldehyde (4% in 0.1 M phosphate buffer). After an overnight post-fixation in paraformaldehyde, the brains were cut in 100 µm thick tangential sections and stained for cytochrome oxidase.

Results

Starting from tangential slices stained for cytochrome oxidase, the reconstructed 2-D barrel map can be obtained in a few clicks by using the provided GUI. The whole procedure takes 7-9 minutes including 3-5 minutes of computation that do not require the intervention of the user. The same process using traditional manual methods (with the help of a raster graphics editor such as Adobe Photoshop) takes, for a well-trained person, 16 to 28 minutes. Most importantly, our automated approach prevents user dependent variability in the obtained barrel map. In order to assess the accuracy of our registration method we first carried out histological control experiments, and then used VSD imaging to functionally evaluate the precision of the barrel map obtained following the full reconstruction procedure.

Histological Validation of the Registration Method

In order to assess the efficiency of our registration method based on automated blood vessel detection and robust ICP, we perpendicularly inserted DiI coated electrodes in the barrel cortex of urethane anesthetized mice (n = 5 experiments, 3 to 6 electrode penetrations per experiments), before processing the brain for histological staining of the cytochrome oxidase following our standard procedures. The electrodes being flat, they did not leave any round or elliptic white marks in the tissue and therefore did not interfere with our registration method (Figure 4). Because the penetration of DiI coated electrodes was perpendicular to the cortical surface, the DiI staining should appear aligned on consecutive cortical sections following proper registration of the slices. In order to control this alignment, the location of DiI spots was reported for each slice on the final fused barrel map image for each section (Figure 4C). The calculated mean distance between DiI spots from the same electrode penetration (34.41 ± 18.93 µm, mean ± SD, n = 5) revealed the subcolumnar resolution of our registration method. These control experiments were further used to evaluate the eventual shrinkage of the cortical tissue due to brain fixation and histological procedures. The distances of 250 µm and 500 µm separating the electrode penetration sites in vivo, were compared to the distances measured between DiI spots on the slices following histological procedures. Over our 5 control experiments, the observed tissue shrinkage within the x-y plane was minimal (< 1.5%).

Assessment of the barrel map reconstruction tool using VSDI of cortical responses to individual whisker deflections

To finally validate our 2-D barrel map reconstruction method, we confronted its resulting map with the functional organization of the barrel cortex established in vivo by real time imaging of cortical responses to individual whisker deflections under urethane anesthesia (n = 4 experiments). Using a mechanical multi-whisker stimulator [START_REF] Jacob | The matrix: A new tool for probing the whisker-to-barrel system with natural stimuli[END_REF], we deflected independently the 24 principal whiskers in a pseudo random order, and imaged the evoked cortical responses in the contralateral barrel cortex using the VSD RH1691. Figure 5 illustrates the results obtained from one experiment. As previously reported in similar conditions [START_REF] Ferezou | Visualizing the Cortical Representation of Whisker Touch: Voltage-Sensitive Dye Imaging in Freely Moving Mice[END_REF], the earliest responses to whisker deflections were localized to the corresponding barrel-related columns (Figure 5A). When reporting the 90% contours of the early cortical responses onto the aligned barrel map, we observe a good anatomo-functional match (Figure 5B). To quantify this match, the distance between the centroid of the anatomically defined barrels and the centroid of the early VSD response (area above a 90% threshold) was measured (Figure 5C). The mean centroid-centroid distance over the 4 control experiments for all the barrels is 60.5 ± 21.2 µm (34.8-116.5 µm range across the barrel field). We observed slightly higher values for the columns located at the border of the map which might result from the curvature of the cortex, the maxima of cortical responses were located within the corresponding barrel area in the majority of cases (86.55%), attesting to the accuracy of the method.

Discussion

We have designed a comprehensive pipeline for the reconstruction of the 2-D barrel map from histological sections. This tool enables a fast reconstruction of a precise barrel map, thus saving a significant amount of time for the experimentalist. Indeed, most of the studies based on optical imaging which required a post-hoc anatomo-functional mapping of the barrel cortex relied on manual reconstruction and alignment of the barrel map with the functional images [START_REF] Kleinfeld | Distributed representation of vibrissa movement in the upper layers of somatosensory cortex revealed with voltagesensitive dyes[END_REF][START_REF] Takashima | Voltage-sensitive dye versus intrinsic signal optical imaging: comparison of optically determined functional maps from rat barrel cortex[END_REF][START_REF] Ferezou | Visualizing the Cortical Representation of Whisker Touch: Voltage-Sensitive Dye Imaging in Freely Moving Mice[END_REF][START_REF] Kerr | Spatial organization of neuronal population responses in layer 2/3 of rat barrel cortex[END_REF][START_REF] Tsytsarev | Study of the cortical representation of whisker directional deflection using voltagesensitive dye optical imaging[END_REF][START_REF] Lustig | Voltage-sensitive dye imaging reveals shifting spatiotemporal spread of whisker-induced activity in rat barrel cortex[END_REF][START_REF] Harris | The effects of focal epileptic activity on regional sensory-evoked neurovascular coupling and postictal modulation of bilateral sensory processing[END_REF]. Other studies [START_REF] Berwick | Integration of neural responses originating from different regions of the cortical somatosensory map[END_REF][START_REF] Berwick | Fine detail of neurovascular coupling revealed by spatiotemporal analysis of the hemodynamic response to single whisker stimulation in rat barrel cortex[END_REF][START_REF] Devonshire | Effects of urethane anaesthesia on sensory processing in the rat barrel cortex revealed by combined optical imaging and electrophysiology[END_REF] used a method based on a warping algorithm described by [START_REF] Zheng | Signal source separation in the analysis of neural activity in brain[END_REF]. However, the manual detection of the fiducial markers required with this approach remains time-consuming in comparison to the automatic detection of blood vessel cross sections proposed here.

In a recent article, Guy and collaborators [START_REF] Guy | Persistence of functional sensory maps in the absence of cortical layers in the somsatosensory cortex of reeler mice[END_REF] used a warping algorithm in order to align the layer IV barrel map, reconstructed manually from histological slices of flattened cortex, with the in vivo functional images using sets of fiducial points. Although the algorithm is not described in detail and probably involves manual selection of the fiducial points, it might be an interesting complement to our approach when working with flattened barrel cortex slices, since it allows a compensation for the curvature of the brain and distortion of the tissue. Finally, instead of using the superficial blood vessels as anatomical landmarks to align the barrel map on the functional images, an alternative approach is to use the images of early cortical responses to single whisker deflections as landmarks [START_REF] Wallace | Plasticity of representational maps in somatosensory cortex observed by in vivo voltage-sensitive dye imaging[END_REF][START_REF] Wang | Voltage-sensitive dye imaging reveals improved topographic activation of cortex in response to manipulation of thalamic microstimulation parameters[END_REF][START_REF] Yang | Thalamic network oscillations synchronize ontogenetic columns in the newborn rat barrel cortex[END_REF]. Using such a method as a standard requires the acquisition of several additional single whisker responses, which might be difficult to implement for instance when working with awake head fixed animals. Note that the tool we propose here to reconstruct the barrel map is valuable whatever solution is chosen in the end to realign the barrel map with the functional images.

On the methodological and mathematical sides, we mainly re-use a set of already existing tools (cross-correlation, ICP and gradient fusion). Our main contribution is to put them together in a coherent processing pipeline. A result of independent interest, that seems to be new, is to show how a family of robust ICP algorithms can be recasted as majorization-minimization algorithms. This in turn allows us to analyze the convergence of these methods.

Obviously the efficacy of the proposed anatomo-functional mapping tool depends upon the quality of the histological slices. Although the preparation of these slices relies on standard protocols which often belong to the daily routines of neurophysiology laboratories, two aspects are essential for the accuracy of the outcome: the quality of the perfusion, and the right thickness of the first (most superficial) slice. Indeed, on the one hand blood vessels have to appear as white circular or elliptical spots on the images to allow the proper registration of consecutive slices and, on the other hand, the superficial blood vessels should be intact to allow the final overlay of the obtained barrel map with the in vivo recordings. When cutting the brain, it is therefore important to set the zero position of the blade with care to ensure a 100 µm thickness to the first slice and thus preserve the integrity of the superficial blood vessels. Finally, although one could deplore that this overlay is a crucial step of the analysis that remains to be achieved manually, we propose here a solution that automatises the most time-consuming phases of the overall process and thus represents a substantial gain in time and precision. Although its use has been demonstrated successfully for VSDI of the adult mouse barrel cortex, it could be expanded to other experimental model species or to the developing brain. Furthermore, the histological section registration method described here might be helpful to reconstruct any tissue in which a majority of blood vessels are orthogonal to the cutting plane of the slices.

Appendix A. Feature Detection with Normalized Cross-Correlation

We consider a set of Gaussian templates {g k } k∈K ,

g k (x) = 1 Z k exp       - ||x|| 2 2σ 2 k      
where Z k is a constant ensuring a normalization x g k (x) 2 = 1. The standard deviations σ k are chosen equally spaced in the range [0, σ max ] (note that σ k = 0 corresponds to a Dirac, supported on a single pixel). Given a template g k , we denote its support as

I k = {x ; g k (x) > η}
where η = 10 -3 is a small tolerance threshold. The normalized cross correlation of a section S with the template g k is then defined as NCC k (S )(x) = y∈I k S (x + y)g k (y)

y∈I k S (x + y) 2 1/2 .
The normalized cross correlation with the whole set of filters is the maximum of all the correlations NCC(S )(x) = max k∈K NCC k (S )(x). (A.1)

A large value of NCC(S )(x) indicates that a vessel is likely to be present at pixel x. In this case, the value k = k(x) of the maximum appearing in (A.1), i.e. such that NCC(S )(x) = NCC k (S )(x), indicates that the radius of this vessel is approximately σ k .

Figure 1 :

 1 Figure1: Cytochrome oxidase staining of tangential sections from the mouse primary somatosensory cortex reveals the structural organization of layer 4 barrels that mirrors the arrangement of the vibrissae on the snout. A. Following the registration of tangential histological slices and the reconstruction of the barrel map, one can see that the spatial organization of the layer 4 barrels matches the layout of the vibrissae on the snout of the animal (B). C. Drawing of a coronal section of the left hemisphere of the mouse brain illustrating the position of layer 4 barrels within the primary somatosensory area of the cortex. After in vivo imaging of barrel cortex activity, sections are cut tangentially to reconstruct the layer 4 barrel map (cutting plane indicated by the red line). D. A series of tangential histological slices stained for cytochrome oxidase. On the first slice one can see superficial blood vessels. On the other slices, one can see white circular to elliptic spots that correspond to sections of plunging blood vessels. Depending on the exact axis of the cut, barrels can be spread over several slices.

Figure 2 :

 2 Figure 2: Alignment of the barrel map with VSDI data using the superficial blood vessels as blueprints.A. Barrel map reconstructed from the histological slices shown in figure1D, using our anatomo-functional mapping method. B. Superimposition of the registered histological slices, allowing the delineation of the superficial blood vessels (in blue). C. The superficial blood vessels also appear in vivo on the fluorescent images taken during the VSDI session. The VSDI data can therefore be aligned with the underlying structural map of the barrel cortex using the blood vessels as anatomical landmarks. On the right, the cortical activity is shown imaged at 10ms following a single C2 whisker deflection with the voltage sensitive dye RH1691 under urethane anesthesia. The barrels outlined from the reconstruction in A are superimposed on the image as white lines.

Figure 3 :

 3 Figure 3: Registration of histological slices and subsequent barrel map reconstruction.A. Snapshots extracted from the GUI illustrating the results of the automated detection of blood vessel cross sections on consecutive histological slices (yellow filled circles and red open circles, respectively -S 3 to S 6 previously shown in Figure1D). The result of the registration obtained by robust ICP applied on the detected blood vessels is illustrated on the right for each pair of slices. Note that the user can delimit manually a region of interest (black dashed lines) in order to restrict blood vessel detection to the barrel field region. B. Images from the same experiment obtained by using our automated barrel map reconstruction tool before (left) or after (right) post-processing steps.

Figure 4 :

 4 Figure 4: Histological validation of the registration method. A. Drawing of a coronal section of the left hemisphere of the mouse brain illustrating the penetration sites of DiI-coated electrodes. B. The DiI fluorescent staining (middle image) does not interfere with the process of blood vessel detection used on bright field images (left), as shown by the absence of detected spots on the DiI stain (right), arrowheads indicate the DiI stained spots corresponding to electrode tracks. C. Following the registration of the slices 2 to 7, we observe a good overlap of DiI spots, indicative of the accuracy of our registration method.

Figure 5 :

 5 Figure 5: Functional validation of the registration method.A. Cortical responses evoked by individual deflections of the 24 principal vibrissae in a urethane anesthetized mouse were imaged using the VSD RH1691. For each whisker, the early cortical response is shown (data averaged from 4 to 18 ms post stimulation time and 40 trials, Gaussian filtered (5 × 5 pxls), and normalized), together with the aligned barrel map obtained from our reconstruction method (white lines). B. The barrel map was extracted from the fused barrel field image (left) and overlaid on the VSD reference image (right) together with the early response locations (90% contours of the early responses (in A) are shown in colors) obtained for each whisker. C. The distance between the centroid of the barrel area and the centroid of the early cortical response (as illustrated in B) was computed, column by column, for the experiment illustrated in A-B (left), and in average for all the control experiments (n = 4).

Table 2

 2 lists the parameters of the method, together with the default values used in our numerical simulations.The successive steps of the algorithms are: Segment the foreground to obtain the input section images {S

	Notation		See
	m ∈ {1, . . . , Q}	Index of a slice	Sec 2.2
	{S 1 , . . . , S Q }	Input slices	Sec 2.2
	NCC	Norm. cross correlation	Appendix A
	X m	Set of detected vessels	Sec 2.2.2
	M	Slices containing barrels	Sec 2.3.2
	T m	Registration map	Eq (5)
	{ S 1 , . . . , S Q }	Registered slices	Eq (5)
	{ S 1 , . . . , S Q }	Inpainted slices	Sec 2.3.1
	S	Fused image	Sec 2.3.2
	S	Drift-corrected result	Eq (6)

1 , . . . , S Q } (Section 2.2.1). For each m ∈ {1, . . . , Q}, compute the list of detected vessel positions X m (Section 2.2.2). For each m ∈ {1, . . . , Q -1}, compute the optimal transform T m to register S m with S m+1 (Section 2.2.3).

Table 1 :

 1 Notations used in the paper.

	Param.		See	Default
	τ	NCC threshold	Sec 2.2.2	0.75
	N	Max. number of vessels	Sec 2.2.2	30
	σ max	Maximum vessel radius	Sec 2.2.2	4.5
	ε	ICP robustness parameter	Eq (2)	-
	σ	Width of debiasing kernel	Eq (6)	10

Table 2 :

 2 Parameters used in the paper.

https://github.com/gpeyre/2014-NeuroMeth-barrels

Acknowledgment

We thank Gérard Sadoc for technical expertise, and Aurélie Daret for her exprimental assistance. This work was supported by Centre National de la Recherche Scientifique (France), the European Research Council (ERC project SIGMA-Vision), the European Union Seventh Framework Programme BrainScaleS (FP7-ICT-2009-6, N 269921), the Agence Nationale pour la Recherche (SensoryProcessing, Transtact).

Appendix B.1. ICP Step 1 with Majorize-Minimize (MM) Iterations

We give here the details of an iterative algorithm to compute a local minimizer (in fact a stationary point) of [START_REF] Egger | 3D reconstruction and standardization of the rat vibrissal cortex for precise registration of single neuron morphology[END_REF], which reads min

This method is similar to re-weighting 2 methods often used for robust ICP (see for instance [START_REF] Bouaziz | Sparse iterative closest point[END_REF]), but we integrate it into a Majorize-Minimize framework, which ensures its convergence.

Starting from an initial transform R (0) , we compute the iterations as

where Ẽ is a so-called surrogate function, which should satisfy (H 1 ) Ẽ(T, T ) -E(T ) is a smooth function of T (of class C 1 ) ;

(H 2 ) for all (T, T ), Ẽ(T, T ) E(T ) ;

(H 3 ) for all T , Ẽ(T, T ) = E(T ).

Under these conditions, it can be shown that the iterations enjoy some good convergence properties. The sequence E(T ( ) ) is decaying and converges to some value E . If E is smooth (which is the case here), ||∇E(T ( ) )|| → 0. Since in our case, the energy E is coercive, the sequence T ( ) is bounded, and all its cluster points T are stationary (i.e. ∇E(T ) = 0) with same energy

The main difficulty in general is to devise a "good" surrogate function Ẽ, i.e. such that one can compute the iteration (B.2) in closed form. The following proposition shows that one can actually design such a surrogate function using a quadratic loss.

2r is decreasing, there exists a constant C(T ) independent of T so that the functional

where

is a majorizing functional for (B.1) and thus satisfies properties

Proof. We rewrite Ẽ as

where we defined ρ(r, r ) = c(r ) + w(r )r 2 where c(r

Thanks to the separability Ẽ (it is a summation over i of functions involving independent variables) and the change of variables r = ||T (x i )z i ||, it thus suffices to prove that ρ is a surrogate functional for ρ on R + . Hypothesis (H 1 ) holds because ρ is C 1 , and one verifies that ρ(r , r ) = ρ(r ) so that (H 3 ) holds.

For any r 0, we consider

It satisfies h(r ) = h (r ) = 0 and h (r) = 2r(w(r )w(r)). Since w is decaying, r is the only point where h is vanishing on R + . This implies that h 0, hence (H 2 ).

The hypothesis that w is decreasing should be interpreted as the condition that ρ should penalize less than a quadratic loss, which makes sense for a robust penalization. Note that the loss (2) that we use in our method satisfies this condition, and that the weighting function satisfies

Appendix B.2. ICP Step 2 with Weighted Quadratic Loss

We consider the problem of solving min

where R is a rotation and t ∈ R 2 . This minimization appears in the MM iteration (B.2) when using the majorizing function (B.3). This problem has a closed form solution, as detailed for instance in [START_REF] Maurer | Registration of 3-d images using weighted geometrical features[END_REF]. For the sake of completeness, we recall the steps of the method. One first centers the points, for i ∈

The optimal rotation is obtained as R = VU T where (U, V) are the eigenvectors of the correlation matrix

(here Λ is the diagonal matrix of eigenvalues). The optimal translation is then computed as

We consider a registered slice S m and its associated vessel locations X m = {x i } i∈I . We denote xi = T m (x i ) the registered vessel locations, where the cumulative transform T m is defined in [START_REF] Land | Cytochrome oxidase staining in the rat smi barrel cortex[END_REF]. We recall that the cross-correlation minimization (as detailed in Appendix A) outputs at each pixel x the index k(x) of the optimal Gaussian template at this location, which has a radius σ k(x) .

We define a mask Φ, which is the set of pixels that are at distance smaller than σ k(x i ) from the point xi . It is thus a union of disks. Pixels in Φ should be discarded and inpainted. This is achieved using a quadratic minimization that seeks a smooth interpolation of missing data where ∇S is a finite difference approximation of the gradient of the image S . The solution of (C.1) corresponds to solving a Poisson equation ∆S = 0 on Φ with Dirichlet boundary conditions given by the constraints. This can be solved using a conjugate gradient method.

Appendix D. Gradient Domain Image Fusion

We consider a set { S m } m∈M of input images to fuse. At each pixel x, we denote the index of largest gradient magnitude as

where ∇ is a finite differences approximation of the gradient operator.

We design a fused vector field as

Since the vector field u is obtained by gluing together gradients from several different images, it is in general not anymore the gradient of an image. We thus reconstruct a valid fused image S using the minimal norm pseudo-inverse, i.e. by computing an image S whose gradient is as close as possible to u with adequate boundary conditions, where ∆ = div •∇ is the Laplacian operator and div = -∇ * is the divergence. When using periodic boundary conditions (which can be used in our case), one solves (D.1) in O(N log(N)) operations (where N is the number of pixels) using an FFT Poisson solver.