Pablo Mesejo
email: pmesejo@ce.unipr.it

Roberto Ugolotti

Ferdinando Di Cunto
email: ferdinando.dicunto@unito.it

Mario Giacobini
email: mario.giacobini@unito.it

Stefano Cagnoni
email: cagnoni@ce.unipr.it

Automatic Hippocampus Localization in Histological Images using Differential Evolution-Based Deformable Models

Keywords: Hippocampus, Deformable Models, Automatic Localization, Histological Images, Global Continuous Optimization, Differential Evolution

In this paper, the localization of structures in biomedical images is considered as a multimodal global continuous optimization problem and solved by means of soft computing techniques. We have developed an automatic method aimed at localizing the hippocampus in histological images, after discoveries indicating the relevance of structural changes of this region as early biomarkers for Alzheimer's disease and epilepsy. The localization is achieved by searching the parameters of an empirically-derived deformable model of the hippocampus which maximize its overlap with the corresponding anatomical structure in histological brain images. The comparison between six real-parameter optimization techniques (Levenberg-Marquardt, Differential Evolution, Simulated Annealing, Genetic Algorithms, Particle Swarm Optimization and Scatter Search) shows that Differential Evolution significantly outperforms the other techniques in this task, providing successful localizations in 90.9% and 93.0% of two test sets of real and synthetic im-

Introduction

Among the different anatomical structures which make up the mammalian brain, the hippocampal formation (HPF) (see Figure 1) is particularly interesting. From an anatomical viewpoint, the HPF, composed by the Hippocampus and the Subiculum (SUB), is located within the medial temporal lobe. In turn, the Hippocampus is composed by the Dentate Gyrus (DG) and Ammon's Horn (CA), which is further composed by three different regions (CA1, CA2, and CA3).

The hippocampus has long been known for its crucial role in learning and memory processes [START_REF] Norman | How hippocampus and cortex contribute to recognition memory: revisiting the complementary learning systems model[END_REF]. Moreover, it has recently been demonstrated that the volume of the hippocampus is an early biomarker for Alzheimer's disease. Therefore, there is a great interest in understanding the cellular and molecular events that take place in this structure, under both normal and abnormal conditions. From this point of view, a precise gene expression map at the cellular and subcellular level within this region can provide crucial information for understanding such biological mechanisms. A very promising data source to derive this map has recently been provided by the Allen Brain Atlas (ABA), a huge, publicly available database that contains highresolution images mapping the expression patterns of most genes contained in 2 the genomes of the analyzed organisms. The first release of the ABA was focused on mouse, and contained the expression patterns of about 20,000 genes obtained by In Situ Hybridization (ISH) of mouse brain sections [START_REF]Allen Reference Atlases[END_REF]. In a more recent release, a similar resource based on spatially mapped microarray data has also been provided for the human brain [START_REF] Hawrylycz | An anatomically comprehensive atlas of the adult human brain transcriptome[END_REF]. There is wide availability of brain images containing morphological and functional information on the hippocampus in different organisms. Thus, it has become extremely important to design image analysis methods that accurately, robustly, and reproducibly identify the hippocampus region, to automatize any relevant analytic procedure. In the ABA, for each gene under consideration, several images are provided, corresponding to different sections of the brain. Each image is labelled according to the corresponding images of a reference atlas comprising 132 coronal and 21 sagittal sections spaced at 100 µm and 200 µm intervals, respectively.

In this paper, we describe a general model-based method for the automatic localization of anatomical districts in histological images, which could be also applied to other structures, as well as to different imaging modalities.

In particular, we describe an application in which such a method is used to locate the hippocampus in ISH images from the ABA. This method allows us to automatically extract image parameters from corresponding regions of a huge number of images, to cluster genes with similar expression patterns and subcellular mRNA distribution. It can be hypothesized that genes whose expression is mapped in images with similar visual features are likely to be In our method, the hippocampus is located by detecting, as landmarks, two regions which are usually well distinguishable within the structure (see Figure 1): the pyramidal (sp) and granule (sg) cell layers, which belong to the CA and DG regions, respectively.

As in most medical imaging applications, the problems we have to solve are mainly related with the characteristics of the images involved in the process. The main problematic features of these images in general, and of the hippocampus region in particular, related with biological features or with the image acquisition process, can be summarized as follows:

• natural variability of brain structure shapes in different subjects;

• fuzziness of the hippocampus boundary;

• limited relevance of color for detecting anatomical structures: regions with similar colors may represent different structures and vice versa, depending on the dye used as well as on local image acquisition settings;

• contrast variability between structures: different genes are not expressed equally in the same anatomical region, making it difficult to construct a consistent model for each landmark in all images. Moreover, grained patterns with many irregularities hamper the classification of individual pixels as belonging to the anatomical structures under consideration;

• orientation issues: the imaged structures may be rotated or displaced on the slice with respect to a "standard" alignment;

• lighting issues: within the same set of images, some are much brighter than others.

• variable resolution even within the same image: high-resolution regions coexist with low-resolution ones;

• presence of artefacts: tears, scraps, bubbles, streaks in tissues, partial cut-off of regions;

• large image size (the typical resolution of ABA images is about 15,000 × 7,000 pixels).

These problems significantly hamper tasks like localization and segmentation of structures in such images.

The fully automatic 2D localization method we propose is based on atlasbased registration and on the optimization of the parameters of a parametric deformable model. Our method can be divided in two stages: (i) selection of the corresponding slice in the reference atlas based on a two-step affine registration, and (ii) proper localization of the hippocampus. This technique was initially presented in [START_REF] Ugolotti | Automatic Hippocampus Localization in Histological Images using PSO-Based Deformable Models[END_REF], where the stochastic search of the hippocampal region was performed using Particle Swarm Optimization (PSO) [START_REF] Kennedy | Particle Swarm Optimization[END_REF]. In this paper, we give more details about the second stage of the method, we delve into the rationale of its operation and we widen its experimental evaluation performing more extensive tests. Moreover, we compare six real-parameter stochastic optimization techniques on the task under consideration; in particular, we show that Differential Evolution (DE) significantly outperforms the other methods taken into consideration: Levenberg-Marquardt (LM) [START_REF] Levenberg | A method for the solution of certain nonlinear problems in least squares[END_REF], Simulated Annealing (SA) [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF], Scatter Search (SS) [START_REF] Glover | Scatter search[END_REF], Genetic Algorithms (GA) [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Ma-chine Learning[END_REF] and PSO. In our tests, we evaluate the localization of the hippocampus in real and synthetic sagittal images, but this method could as well be applied to other subcortical structures, image modalities or anatomical planes.

The paper is organized as follows: in section 2 we provide the theoretical foundations of our work, as well as an overview of previous related work. In section 3, a general overview of the method is presented, providing details about the implementation of the localization of the hippocampus. Finally, section 4 presents results on both synthetic and real images with appropriate statistical tests, followed, in section 5, by some final remarks and a discussion about possible future developments. Figure 2: Hippocampus variability. Horizontal-wise: hippocampus sections taken at different levels look different, and different genes can as well produce very different visual features. Vertical-wise: corresponding sections from different brains maintain some shape similarity.

Theoretical background and related work

Deformable Models

Deformable models are curves or surfaces, defined within the image domain, that are deformed under the influence of "internal" forces, related with the curve features, and "external" forces, related with the surrounding image. Internal forces enforce regularity constraints and keep the model smooth during deformation, while external forces are defined such that the model is attracted toward an object or other features of interest within the image.

The term "deformable models" was first used in the late eighties [START_REF] Terzopoulos | Deformable Models[END_REF][START_REF] Terzopoulos | Constraints on Deformable Models: Recovering 3D Shape and Nonrigid Motion[END_REF], and one of the first examples, called "snakes" or Active Contour Models, was presented shortly after in [START_REF] Kass | Snakes: Active contour models[END_REF].

Active Shape Models (ASMs) [START_REF] Cootes | Active shape models-their training and application[END_REF] add more prior knowledge to deformable models. These shape models derive a "point distribution model" from sets of labelled points (landmarks) selected by an expert in a training set of images.

Thus, the final model of the object of interest can be derived by examining the distribution of the positions of the labelled points. The model considers the points' average positions and the main modes of variation observed in the training set. It is important to notice that an instance of the model can only take into account deformations which appear in the training set: because of this, the model may have problems with unexpected shapes, but it is robust with respect to noise and image artefacts, like missing or damaged parts.

Although originally developed for computer vision applications to natural scenes and computer graphics problems, the applicability of deformable models in medical image analysis has already been proven [START_REF] He | A comparative study of deformable contour methods on medical image segmentation[END_REF][START_REF] Heimann | Statistical shape models for 3D medical image segmentation: a review[END_REF].

Differential Evolution

In the last two decades, research on global optimization has been very active [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF][START_REF] Herrera | Tackling Real-Coded Genetic Algorithms: Operators and Tools for Behavioural Analysis[END_REF][START_REF] Deb | A computationally efficient evolutionary algorithm for real-parameter optimization[END_REF], and many different deterministic and stochastic algorithms for continuous optimization have been developed. Among the stochastic approaches, Evolutionary Algorithms (EAs) [START_REF] Eiben | Introduction to Evolutionary Computing[END_REF][START_REF] Engelbrecht | Computational Intelligence: An Introduction[END_REF] offer a number of advantages that make them attractive: implicit parallelism, robust and reliable performance, global search capability, no need of specific information about the problem to solve, easy implementation, good insensitivity to noise, and no requirements for differentiable or continuous objective functions. DE, first introduced by Storn and Price [START_REF] Storn | Differential Evolution-A Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces[END_REF], has recently been one of the most successful evolutionary algorithms. Unlike traditional EAs, DE perturbs the current population members with the scaled differences of randomly selected and distinct individuals [START_REF] Das | Differential Evolution: A Survey of the Stateof-the-Art[END_REF]. This way, in the first iterations, the elements are widely scattered in the search space and have a great exploration ability.

During optimization, the individuals tend to concentrate in the regions of the search space with the best fitness values [START_REF] Neri | Recent advances in differential evolution: a survey and experimental analysis[END_REF].

In DE, new individuals that will be part of the next generation are created by combining members of the current population. Every individual acts as a parent and is associated to a donor vector. In the basic version of DE, the donor vector V i for the i th parent (X i ) is generated, by combining three random and distinct population members (X r1 , X r2 ) and X r3 , as follows:

V i = X r1 + F • (X r2 -X r3 )
where F (scale factor) is a real-valued parameter that strongly influences DE's performances and typically lies in the interval [0.4, 1]. Other mutation strategies have been applied to DE, experimenting with different base vectors and different numbers of vectors for perturbation.

After mutation, every parent-donor pair generates an offspring O i by means of a crossover operation. The newly-generated offspring is evaluated and its fitness and its parent's are compared. The better survives and will be part of the next generation.

Related Work

We shortly review some of the most relevant and recent approaches which combine deformable models and evolutionary computation techniques. In [START_REF] Ballerini | Genetic snakes: active contour models by genetic algo-rithms[END_REF], "Genetic snakes" are active contour models with an optimization procedure based on genetic algorithms. In [START_REF] Mcintosh | Medial-based Deformable Models in Nonconvex Shape-spaces for Medical Image Segmentation using Genetic Algorithms[END_REF] a GA evolves a population of medial-based shapes, using prior shape knowledge to produce feasible deformations while also controlling the scale and localization of these deformations. In [START_REF] Ghosh | A Genetic Algorithm-Based Level Set Curve Evolution for Prostate Segmentation on Pelvic CT and MRI Images[END_REF] a GA is used to perform level set curve evolution using texture and shape information to automatically segment the prostate in CT and MRI pelvic images. Finally, in [START_REF] Ibáñez | Genetic approaches for topological active nets optimization[END_REF], the authors describe a GA-based method that minimizes the Topological Active Nets energy, evaluating the segmentation results of a greedy, a genetic, and a memetic algorithm over seven synthetic images and three CT scans of human bones.

Recently, PSO has been successfully used in conjunction with various types of deformable models. Asl and Seyedin [START_REF] Asl | Active Contour Optimization using Particle Swarm Optimizer[END_REF] apply the technique proposed in [START_REF] Maceachern | Genetic algorithms for active contour optimization[END_REF] using PSO instead of a GA, obtaining similar results in terms of precision but in shorter time. In [START_REF] Feltell | 3D level set image segmentation refined by intelligent agent swarm[END_REF] the initial segmentation based on the level set method is refined using swarms of intelligent agents. Finally, in [START_REF] Shahamatnia | Application of particle swarm optimization and snake model hybrid on medical imaging[END_REF], a customized PSO algorithm overcomes the drawbacks of snakes, including initialization, concave boundaries, sensitivity to noise and local minima.

To the best of our knowledge, one of the few works that uses DE and deformable models for solving medical imaging problems is presented in [START_REF] Novo | Topological Active Models optimization with Differential Evolution[END_REF],

where DE, in combination with a greedy search algorithm, is used to evolve Topological Active Nets, a discrete implementation of an elastic mesh, for CT image segmentation.

DE-based hippocampus localization

Our accurate and automatic structure localization method consists of two phases:

1. Best Reference Slice Selection, implemented as a two-step affine registration method which: (i) determines the position of the section displayed in the target image according to a reference atlas, and (ii) extracts the region of interest (ROI) where the hippocampus is more likely to be located. Hence, it represents the initialization of the deformable model: it determines which template should be used and where it should be applied.

Best Reference Slice Selection

The initialization of the Deformable Model, in terms of choosing a model and its starting position, is obtained using an atlas-based affine registration with the reference images of the ABA. The main idea is to find the sagittal reference slice of the atlas which best matches the target image. This phase produces two results: firstly, based on the information contained in the corresponding reference atlas image, it allows one to extract the Region of Interest (ROI) where the hippocampus is expected to be located. Secondly, it makes it possible to determine the position, within the brain, of the section represented in the target image. Consequently, one can select the corresponding point distribution model of the hippocampus (derived empirically as described for the ASMs), which is to be applied in the following step within the selected ROI. This can be achieved by the two-step matching method presented in [START_REF] Ugolotti | Automatic Hippocampus Localization in Histological Images using PSO-Based Deformable Models[END_REF].

The first step performs a "global" affine registration between the target image and each image in the reference atlas, considering the shapes of the whole sections and ignoring the internal structures. Then, using the atlas, the ROI where the hippocampus is more likely to be located is extracted, and a "local" affine registration is performed in order to match the hippocampus in the reference and the target image. After these two registration steps, two similarity measures between the target image and each reference image are computed. The reference image that achieves a better trade off between these two measures is selected as the best reference.

In order to evaluate the correctness of this two-step approach, we computed the distance (in slices) between the slice selected by our system and the one suggested by the ABA on a test set of 320 images. In 45% of cases the slice selected was exactly the same, while in 43%, 10% and 2% of cases there was a difference of one, two, and three or more slices, respectively. In those cases in which our system selected a reference slice different from the one proposed by the ABA, we also compared the fitness values obtained by the hippocampus localization algorithm. The results of these experiments showed that there are no statistically significant differences between the two choices. Therefore, these results demonstrate that our initialization method is effective and can be used also in problems where no ground truth is available.

Hippocampus Localization

While in the standard ASMs external forces are driven by the contours in the image, in our approach we parametrically deform the model shape to match as closely as possible the shape of the hippocampus in the region we want to locate. The model is moved and deformed, by altering its parametric representation using an optimization heuristic like DE or PSO, which maximizes a function which measures the similarity between the model and the object itself. Since, in our problem, color and shapes are greatly variable, the only common information on which we can rely is that all hippocampi have two substructures (sg and sp) whose shapes have more consistent features, besides being usually characterized by lower intensity values and higher color saturation with respect to the surrounding structures.

To enhance the contrast between sg, sp and the surrounding parts, we pre-process the ROI containing the hippocampus by stretching its histogram.

We invert pixel intensities, to consistently use the convention by which, as for binarized images, the most relevant regions in the image have higher intensities. These operations are not applied to the original images of the ABA we show, to avoid altering their actual appearance.

Templates

Obviously, our main concern is to adapt a model that permits an efficient localization, but we want, as well, to obtain this goal with the simplest possible model, i.e. having the least possible number of points in each template. In this work, for every slice of the reference atlas, two templates (one composed of 7 points for sg, one with 8 points for sp) have been created by manually selecting significant corresponding landmarks in all images in a training set.

A template aims at taking into account all possible positions and deformations with respect to the prototypical shape of the structure it represents.

The template does not refer to the absolute positions of the points, but describes the relative positions (or shifts) between consecutive points in polar coordinates. This is uncommon for a DM, since the majority of the implementations use cartesian coordinates. Every template is composed of two parts:

• an "inner set" of points that lies on the anatomical part we want to locate;

• an "outer set" of points that lies just outside it, obtained by rigidly shifting the previous set.

For a proper localization, the "outer set" is fundamental; otherwise, using only the "inner set", a completely dark and large area would always achieve the highest values of the target function.

A template is fully described by four 2×n matrices, where n is the number of points in the template:

S =          ρ m1 ϑ m1 ρ m2 ϑ m2 . . . . . . ρ mn ϑ mn          ∆ =          ∆ ρ1 ∆ ϑ1 ∆ ρ2 ∆ ϑ2 . . . . . . ∆ ρn ∆ ϑn          L =          ρ l1 ϑ l1 ρ l2 ϑ l2 . . . . . . ρ ln ϑ ln          U =          ρ u1 ϑ u1 ρ u2 ϑ u2 . . . . . . ρ un ϑ un         
S is the "standard template" and represents the standard coordinates of the inner set, ∆ is the displacement of the outer set with respect to the inner set, L and U are the minimum and maximum values allowed for every parameter that describes the inner set. It should be noticed that ρ 1 and ϑ 1 represent the positions of the first point with respect to the upper left corner of the image. After that point, proceeding row-wise, every (ρ, ϑ)

pair represents the shifts of the subsequent point in terms of distance and angle with respect to the previous point, respectively. As a consequence, the number of parameters to optimize for each model is twice the number of the points in the template (14 parameters for sg and 16 for sp).

The matrices S, L and U have been computed by manually selecting the reference points in a training set. The first one is the median of the selected shifts and the other two are the minimum and the maximum values observed in the training set, respectively. To improve the templates, a manual refinement of the parameters has been performed. The matrix ∆ has been manually built based on the observation of several hippocampi.

From a computational point of view, a model is a vector of 2 × n elements that ranges within the "boundary templates" L and U; its elements (coordinates) are evolved by a metaheuristic which optimizes the target function described in section 3.3. The aims with which the target function was defined are to match the model with the target while keeping it as close as possible to the shape of the standard template S.

Target Function

As in classical deformable models, we consider the model to be subjected to external forces (driven by the image features) and internal forces (driven by the model itself).

The target function F , which is to be maximized, has three components:

external energy E, internal energy I, and contraction factor C:

F = E -(I + C)
The external forces move (and deform) the model to maximize the intensity of pixels in the inner set, while minimizing the intensity of pixels in the outer set. For both sets, we evaluate the intensity of the image within a 3×3 neighborhood N 3 of all points in the model (Punctual Energy, P E) and in p intermediate points along the segment between two consecutive points (Continuous Energy, CE), i.e.

P E = n i=1 [T (N 3 (I i )) -T (N 3 (O i ))]
where n is the number of points in the model,

I i = {x i , y i }
T (I i-1 + j p + 1 (I i -I i-1 )) - n i=2 p j=1 T (O i-1 + j p + 1 (O i -O i-1 ))
where p is the number of points to evaluate in every segment. In our case we set p = 20. The final external energy is computed as

E = γ P • P E + γ C • CE
where the weights γ P and γ C have been empirically set to 5 and 1 for all tests, respectively (see Figure 3).

The internal energy I, related with the forces that oppose the deformation of the model, is computed as:

I = ξ ρ • n i=2 (ρ i -ρ mi ) 2 + ξ ϑ • n i=2 (ϑ i -ϑ mi ) 2
where ξ ρ and ξ ϑ are two positive weights that regulate the deformability of the model. The higher their values, the less the model deforms and the more it keeps the shape of the standard template S. Please notice that the index i starts from 2. This way, the first shift (which actually represents the starting position of the model) is not taken into account and the model is independent of its position in the image.

Finally, the contraction factor C also regulates the model's deformability to avoid unfeasible situations (e.g. having the extremes of sg too close to each other), and is defined as follows:

used in [START_REF] Ugolotti | Automatic Hippocampus Localization in Histological Images using PSO-Based Deformable Models[END_REF], where it was shown to give better results than the standard PSO on this task; SS with local search based on the line search algorithm and the blx-α crossover, as described in [START_REF] Herrera | Continuous scatter search: An analysis of the integration of some combination methods and improvement strategies[END_REF]; GA and SA. In LM and SA we used, as initial value, the best solution in a randomly generated set of the same size as the other methods' populations.

The parameters (as well as the population topologies and the mutation and crossover operators) set in the tests have been initially chosen based on the most commonly used, and then refined during the development of the system. The values chosen for the most relevant parameters are shown in Table 1. For all techniques, a limit of 20000 function evaluations has been set as the only termination criterion, and the limits for the deformation were used as constraints in creating new solutions by crossover and mutation.

22

Table 2 summarizes the results of our tests. In the upper part of the table we report the results for the 20 synthetic images and, in the lower part, the ones related to the 320 real images. The second column reports the average fitness for all tests, while the third reports the average standard deviation of fitness in the 25 experiments performed on each image. The fourth and fifth columns report the average of the worst and best results obtained for each image. Similarly, the sixth column reports the average of the medians.

We focused our analysis on the comparison between the performance of DE and of the other methods. DE achieved the best average fitness as well as the lowest standard deviation, which indicates a more robust behavior, as it is able to produce more consistent results over different trials; the performance of the other algorithms is more dependent on their random initialization. To assess the significance of this result, we performed a statistical test with confidence level of 0.01. We used non-parametric tests because the assumption of normality was not met. The Kruskal-Wallis test [START_REF] Kruskal | Use of Ranks in One-Criterion Variance Analysis[END_REF] was statistically significant in all cases (p-value < 1.00E-16) and proved the existence of differences between the sets of results, where at least one sample median was significantly different from the others. After that, the paired Wilcoxon signed-rank test was performed [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF], assuming as null hypothesis that the median of DE results is less than or equal to the median of the other methods. This is a one-tailed test in which the alternate hypothesis is that the median of any other method is less than the median of DE. The corre-sponding p-values (shown in the last column of Table 2) reject, in all cases, the null hypothesis, showing that significant differences exist between the performance of DE and of the other methods. Moreover, the average worst values show that DE avoids local minima and is able to get good results more consistently than the other methods. than the other techniques but, while the former continues to improve its performance, the latter is not able to refine the solution found. At the same time, the evolution of SS is more discontinuous due to the premature convergence of the reference set, followed, after each reconstruction, by an abrupt improvement of the fitness values. The worst results were obtained by the gradient-based local search method (LM), due to the multimodality and non-differentiability of the function to optimize, providing a reason for using metaheuristics to solve this kind of problems.

Overall System Performance

In a second test we evaluated the results of the entire localization method by dividing the outcomes into three quality classes (see Figure 5): which must be detected, but they do not cover them entirely or (ii) at most two points are slightly outside of them;

3. Error: all other possibilities, from three or more misplaced points to models which are located in a completely different position of the brain. If the results obtained by the model chosen in the Best Reference Slice Selection phase were not good enough (using a threshold for the target function values), we repeated the procedure using the model which ranked second in the previous phase, and took the best result.

Our method was able to perfectly localize the hippocampus in 58.0% of synthetic images and in 47.8% of real images, and reached a good localization in 35% and 43.1% of cases, respectively. This means that our method was able to localize the hippocampus satisfactorily in 93.0% with synthetic images (20 images and 25 runs per image) and in 90.9% of cases with real images (320 images and 1 run per image).

The most problematic situations that appear to drive the system to bad localizations are substantially three:

• low quality of images (for instance, images that are damaged or have very low contrast): this can affect the results of ROI extraction or of the Best Reference Slice Selection phase;

• images where the hippocampus has low levels of gene expression and is surrounded by other anatomical parts with higher expression levels;

• hippocampi having shapes which differ substantially from the typical cases included in the training set.

Most of these errors could most probably be avoided by increasing the size of the training set for the point distribution models, or improving the preprocessing phase of the images.

On the other hand, the method we adopted shows very good performances when dealing with corrupted images of hippocampi.

With regard to time complexity, the localization task, implemented in C, employs on average 2.88 seconds per image, applying the Best Reference Slice Selection process on 10 reference slices.

Discussions and Conclusions

In this work we presented a two-step algorithm aimed at automatically localizing the hippocampus in histological images. The first phase roughly locates the anatomical structure by comparing the image under consideration with images taken from the ABA reference atlas. In the second phase, the structure is more precisely localized using a DE-based parametric deformable model that adapts its shape to match the anatomical structure of interest.

The method is able to deal with imprecise and incomplete images, and, in our tests on actual images from the ABA, has been successful in 90.9% of cases.

Since classical gradient-based local search methods (like LM) are not able to solve this problem satisfactorily, due to its multimodality and high dimensionality, we have applied and studied different metaheuristics: DE, SA, GA, PSO, and SS. DE achieved the best results, for both average and standard deviation. This result can be explained by a better balance between exploration/diversification and exploitation/intensification and upholds the results obtained in other cases by DE on benchmarks containing multimodal nonseparable functions [START_REF] Vesterstrom | A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms[END_REF]. Moreover, DE has shown greater robustness: it is able to reproduce good results more consistently over many trials, whereas the performance of other algorithms, like GA or SS, is more dependent on the stochastic initialization of individuals and parameters. DE is also, together with GA, the metaheuristic that converges faster to good solutions.

Possible future work includes automating the configuration step which sets the algorithm parameters. Another interesting development could consider applying this method to other subcortical structures and anatomical planes in order to evaluate its general applicability. The automatic segmentation of the structure and the texture analysis of the segmented areas will lead, in a subsequent phase of the project within which this method was developed, to identifying sets of genes whose expression generates similar textural patterns in corresponding regions. We aim to cluster similar genes according to such visual features to detect correlations in their expression, since it can be argued that genes with similar expression patterns are also functionally correlated.

putation, CEC '04, 1980-1987, 2004. 605

Figure 1 :

 1 Figure 1: Regions in Hippocampal Formation (HPF)

2 .

 2 Structure Localization, in which the point distribution model selected in the previous step is adapted by a stochastic optimization procedure based on DE, to fit the image region where the structure of interest (if any) is located.

Figure 3 :

 3 Figure 3: Inner and outer models superimposed to a hippocampus image, showing the components of the external energy.

  schedule: T k+1 = 0.8•T k DE/target-to-best/1 Population Size = 80 Tournament of size 4 Population Size = 80 Iterations = 250 Uniform Crossover Iterations = 250For every image, we first performed the extraction of the ROI and the selection of the best model, as described in Section 3.1. After this, for each optimization technique, we ran 25 tests on every image, for a total of 8000 experiments on real images and 500 on synthetic images.

Figure 4

 4 Figure 4 describes the behavior of the five metaheuristics and a gradientbased local search method by plotting the average fitness value versus the number of function evaluations. It can be seen how DE and GA start better

Figure 4 :

 4 Figure 4: Evolution of the fitness function for the localization of sp (above) and sg (below).

Figure 5 :

 5 Figure 5: Upper row: perfect matches; middle row: good results; lower row: erroneous localizations.

  is the i th point of the inner set (in cartesian coordinates), O i = {x i + ∆x i , y i + ∆y i } is the i th point of the outer set, T (P ) is the intensity of the image in P , if P is a point, or the average intensity, if P is a neighborhood; ∆x i , ∆y i are the elements of ∆ in cartesian coordinates. As well,

	n	p
	CE =
	i=2	j=1

Table 1 :

 1 Parameters used in testing different optimization techniques.

Table 2 :

 2 Comparative results of sp and sg localization for synthetic and real images (20000 function evaluations per experiment; the higher the values, the better the localization).

				Synthetic Images
				sp localization
	Method Average StdDev Avg. Worst Avg. Best Avg. Median Wilcoxon test
	DE	120.41	5.49	106.30	126.37	121.68	-
	PSO	105.86	8.67	84.86	118.61	107.71	<1.00E-16
	SS	105.70	6.16	91.66	115.69	106.37	<1.00E-16
	GA	105.93	15.64	54.93	120.40	110.11	<1.00E-16
	SA	101.80	8.19	80.87	110.64	103.76	<1.00E-16
	LM	35.43	30.86	19.96	94.24	34.11	<1.00E-16
				sg localization
	Method Average StdDev Avg. Worst Avg. Best Avg. Median Wilcoxon test
	DE	128.91	2.92	121.18	132.99	129.49	-
	PSO	118.29	5.85	105.55	127.63	119.05	<1.00E-16
	SS	115.56	6.00	101.37	124.81	116.08	<1.00E-16
	GA	115.03	9.67	87.86	127.20	116.49	<1.00E-16
	SA	111.23	3.87	102.49	117.91	111.40	<1.00E-16
	LM	64.76	19.58	33.79	103.34	61.42	<1.00E-16
				Real Images
				sp localization
	Method Average StdDev Avg. Worst Avg. Best Avg. Median Wilcoxon test
	DE	141.47	5.59	124.36	146.06	145.21	-
	PSO	133.56	8.24	110.94	143.89	138.64	<1.00E-16
	SS	131.06	6.04	118.54	141.66	135.81	<1.00E-16
	GA	132.00	14.49	86.07	143.70	138.58	<1.00E-16
	SA	131.00	3.05	123.83	137.71	134.91	<1.00E-16
	LM	79.32	28.26	9.60	126.11	85.02	<1.00E-16
				sg localization
	Method Average StdDev Avg. Worst Avg. Best Avg. Median Wilcoxon test
	DE	148.24	1.54	143.58	149.73	149.30	-
	PSO	144.65	3.69	135.37	148.70	145.52	<1.00E-16
	SS	140.07	4.27	129.10	146.23	140.69	<1.00E-16
	GA	140.89	6.83	122.56	148.50	143.06	<1.00E-16
	SA	138.31	1.70	134.60	141.51	139.28	<1.00E-16
	LM	105.58	20.01	65.54	141.61	110.52	<1.00E-16

Acknowledgments

Pablo Mesejo is funded by the European Commission (MIBISOC Marie Curie Initial Training Network, FP7 PEOPLE-ITN-2008, GA n. 238819). This work was also funded by Compagnia di San Paolo, Torino, Italy (Neuroscience Programme). All the mouse brain images were downloaded from the website: Allen Mouse Brain Atlas [http://mouse.brain-map.org]. Seattle (WA): Allen Institute for Brain Science. c 2009.

If ξ c < 0 the two extremes of the model repel each other, if ξ c > 0 they attract each other. In our case we set ξ c > 0 for the sp models and ξ c < 0 for the sg models.

Experimental Results

To test the system we used both real and synthetic images. As for the former, we randomly selected 320 images (corresponding to 320 different genes) from the ABA within subsets of representative samples of all possible hippocampi, featuring high and low levels of gene expression, good-quality and low-quality images, and so on. We also created 20 synthetic images that represent simplified versions of the real ones. In these images, representing all reference slices, the hippocampus is made up of small circles having random radius and color; small and big ellipses were also added trying to simulate cells and, finally, gaussian (µ ∈ [0.0, 0.2] and σ ∈ [0.01, 0.15]) and salt and pepper noise (density ∈ [0.05, 0.25]) were introduced to add fuzziness to the images.

Tests were run on a computer equipped with a 64-bit Intel R Core TM i7 CPU running at 2.67 GHz with 4 Gb of RAM.

We compared the results obtained by DE with those obtained by LM, a classical combination of the gradient descent and the Gauss-Newton methods, and the following stochastic optimization techniques: a variation of PSO,


Roberto Ugolotti

Pablo Mesejo

Samantha Zongaro

Barbara Bardoni

Gaia Berto

Federico Bianchi

Ivan Molineris

Mario Giacobini

Stefano Cagnoni
email: cagnoni@ce.unipr.it

Ferdinando Di Cunto
email: ferdinando.dicunto@unito.it

Visual Search of Neuropil-Enriched RNAs from Brain In Situ Hybridization Data through the Image Analysis Pipeline Hippo-ATESC

Motivation: RNA molecules specifically enriched in the neuropil of neuronal cells and in particular in dendritic spines are of great interest for neurobiology in virtue of their involvement in synaptic structure and plasticity. The systematic recognition of such molecules is therefore a very important task. High resolution images of RNA in situ hybridization experiments contained in the Allen Brain Atlas (ABA) represent a very rich resource to identify them and have been so far exploited for this task through human-expert analysis. However, software tools that may automatically address the same objective are not very well developed.

Results: In this study we describe an automatic method for exploring in situ hybridization data and discover neuropilenriched RNAs in the mouse hippocampus. We called it Hippo-ATESC (Automatic Texture Extraction from the Hippocampal region using Soft Computing). Bioinformatic validation showed that the Hippo-ATESC is very efficient in the recognition of RNAs which are manually identified by expert curators as neuropil-enriched on the same image series. Moreover, we show that our method can also highlight genes revealed by microdissection-based methods but missed by human visual inspection. We experimentally validated our approach by identifying a non-coding transcript enriched in mouse synaptosomes. The code is freely available on the web at http://ibislab.ce.unipr.it/software/hippo/.

Introduction

The communication between neuronal cells is primarily achieved through chemical transmission at synapses, specialized subcellular structures in which axons and dendrites of connected neurons are closely juxtaposed. In the case of excitatory contacts, the most intensely studied and better understood type, synapses are formed at structures known as dendritic spines (DS), small protrusions of the dendritic membrane that compartmentalize the biochemical events activated by synaptic transmission [START_REF] Hering | Dendritic spines: structure, dynamics and regulation[END_REF][START_REF] Tada | Molecular mechanisms of dendritic spine morphogenesis[END_REF]. One of the most remarkable features of DS is that their shape and efficiency can be individually modified as a function of activity [START_REF] Saneyoshi | Regulation of spine and synapse formation by activity-dependent intracellular signaling pathways[END_REF]. Targeting of coding and non-coding RNAs to axons, dendrites and to dendritic spines (collectively referred to as the neuropil) plays a very important role in the localized control of gene expression underlying these phenomena [START_REF] Mikl | RNA localization in neurite morphogenesis and synaptic regulation: current evidence and novel approaches[END_REF][START_REF] Swanger | Making and breaking synapses through local mRNA regulation[END_REF][START_REF] Tolino | RNA-binding proteins involved in RNA localization and their implications in neuronal diseases[END_REF]. Among the protein-coding RNAs, Ca2+-calmodulin-dependent protein kinase alpha subunit (Camk2a) [START_REF] Burgin | In situ hybridization histochemistry of Ca2+/calmodulin-dependent protein kinase in developing rat brain[END_REF][START_REF] Mayford | The 39-untranslated region of CaMKII alpha is a cis-acting signal for the localization and translation of mRNA in dendrites[END_REF], Map2 [START_REF] Garner | Selective localization of messenger RNA for cytoskeletal protein MAP2 in dendrites[END_REF], Shank [START_REF] Bockers | Differential expression and dendritic transcript localization of Shank family members: identification of a dendritic targeting element in the 39 untranslated region of Shank1 mRNA[END_REF], b-actin [START_REF] Tiruchinapalli | Activity-dependent trafficking and dynamic localization of zipcode binding protein 1 and beta-actin mRNA in dendrites and spines of hippocampal neurons[END_REF] and Arc [START_REF] Bramham | Local protein synthesis, actin dynamics, and LTP consolidation[END_REF] are the best documented examples of neuropilenriched mRNAs. On the other hand, the dendrite enriched noncoding transcript Bc1 has been shown to regulate synaptic plasticity by locally repressing the translational of mGluR receptors [START_REF] Zhong | BC1 regulation of metabotropic glutamate receptor-mediated neuronal excitability[END_REF]. Therefore, the automatic detection of neuropilenriched transcripts has become a very important issue. This problem has been so far addressed through two main general strategies: the direct measurement of RNA molecules in extracts prepared from microdissected neuropil or the use of RNA in situ hybridization on cells and tissues. The rodent hippocampus has been fruitfully used for microdissection-based expression studies, in virtue of the very precise arrangement of cell bodies and neuronal projection that characterize this structure [START_REF] Poon | Identification of process-localized mRNAs from cultured rodent hippocampal neurons[END_REF][START_REF] Zhong | Dendritic mRNAs encode diversified functionalities in hippocampal pyramidal neurons[END_REF]. In combination with microarray-based measurements, this technique has led to the identification of approximately 200 transcripts enriched in the neuropil, as compared to cell bodies [START_REF] Poon | Identification of process-localized mRNAs from cultured rodent hippocampal neurons[END_REF][START_REF] Zhong | Dendritic mRNAs encode diversified functionalities in hippocampal pyramidal neurons[END_REF][START_REF] Cajigas | The local transcriptome in the synaptic neuropil revealed by deep sequencing and highresolution imaging[END_REF]. More recently, an RNA-seq-based study has shown that the hippocampal neuropil contains approximately 2550 coding transcripts [START_REF] Cajigas | The local transcriptome in the synaptic neuropil revealed by deep sequencing and highresolution imaging[END_REF]. However, such a study did not provide information about the relative enrichment of transcripts between the neuropil and the cell bodies and did not report any information about the non-coding transcripts.

Systematic in situ hybridizations on adult mouse brain with probes derived from virtually all protein coding genes and from many non-coding transcripts have been performed within the Allen Brain Atlas project (ABA) [START_REF] Lein | Genomewide atlas of gene expression in the adult mouse brain[END_REF][START_REF] Ng | Large-scale neuroinformatics for in situ hybridization data in the mouse brain[END_REF]. The analyses of ABA images so far performed have identified many neuropil-enriched transcripts [START_REF] Lein | Genomewide atlas of gene expression in the adult mouse brain[END_REF], but suffer of two main limitations. On the one hand, although the resolution limit of the ABA images is 1.07 mm, theoretically allowing for the discrimination of sub-cellular details, their systematic mining by automatic tools has been focused on exploring the general gene expression patterns at low resolution (200 mm), or on measuring expression in cell bodies [START_REF] Lein | Genomewide atlas of gene expression in the adult mouse brain[END_REF]. On the other hand, human expert inspection of high-resolution images has led to the highly specific identification of 57 dendrite-enriched RNAs, but may have significantly underestimated the number of neuropil-enriched RNAs, as it would seem if considering the much higher number of neuropil transcripts detected by RNA-seq [START_REF] Cajigas | The local transcriptome in the synaptic neuropil revealed by deep sequencing and highresolution imaging[END_REF].

In this study we describe the implementation of an automatic pipeline aiming at detecting transcripts enriched in the hippocampal neuropil of adult mice, by systematically exploring the high resolution images contained in the ABA. The method is based on the automatic identification of the different hippocampal sub-regions in high resolution ABA images, followed by the extraction of analysis of many different image-texture features. On this basis, we ranked the mouse coding and non-coding transcripts represented within the ABA according to their similarity to well known neuropil-enriched transcripts. The comparison of our ranking with the results of microdissection studies confirmed the high specificity of our method. We experimentally validated our results by identifying a new non-coding transcript associated to the synaptodendrosomal compartment.

Methods

The Hippo-ATESC Pipeline

The automatic pipeline is based on three main steps: i) localization of relevant hippocampal sub-regions; ii) characterization of the texture of these regions; iii) training of a model for neuropil-enriched transcripts, on the basis of prototype mRNAs. A schematic representation of the procedure is given in Fig. 1. In order to identify the main hippocampal regions within ABA in situ hybridization images, we adapted an automatic method that we previously described [START_REF] Mesejo | Automatic hippocampus localization in histological images using Differential Evolutionbased deformable models[END_REF]. In this method, the localization of hippocampal main structures, i.e. the Ammon's Horne (AH) and the dentate gyrus (DG) is achieved by searching the parameters of an empirically-derived deformable model (DM) [START_REF] Jain | Deformable template models: a review[END_REF], which maximizes its overlap with the corresponding anatomical structure in the brain image using a metaheuristic [START_REF] Gendreau | Handbook of Metaheuristics[END_REF] and gives as output the location of the hippocampal region achieving a good trade-off between accuracy and speed [START_REF] Mesejo | Automatic hippocampus localization in histological images using Differential Evolutionbased deformable models[END_REF] (Fig. 2A). Afterwards, we performed a segmentation step in order to identify 13 different hippocampal sub-regions, corresponding to specific sub-sections of the dentate gyrus and of the CA1 and CA3 regions of the Ammon's Horn (Fig. 2B). From these regions, we extracted 220 textural features from each image. From this large dataset we first extracted a training set composed of three prototypical neuropilenriched genes and 17 negative examples, based on which a genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] selected a subset of 52 significant features (Table S1). Finally, for every gene, we extracted its 52-feature vector and compared it to a reference vector (the average of the three prototype vectors) using the Pearson coefficient as a distance measure [START_REF] Rodgers | Thirteen Ways to Look at the Correlation Coefficient[END_REF].

Localization and Segmentation of Hippocampal AH and DG

The method used to localize the AH and the DG of the hippocampus within ABA sagittal sections is described in more details in [START_REF] Mesejo | Automatic hippocampus localization in histological images using Differential Evolutionbased deformable models[END_REF]. Briefly, it can be divided into two stages: (i) selection of the corresponding slice in the reference atlas based on a twostep affine registration process, and (ii) proper localization of the hippocampus. The first stage recognizes the location, within the brain, of the section displayed in the target image based on a comparison with the images of a reference atlas, and extracts the region of interest (ROI) where the hippocampus is most likely to be located. The following stage performs the proper localization of the anatomical structure under study using statistical shape models. In particular, before localization, the ROI that has been extracted is preprocessed to remove noise and refine the detection of the background. Consequently, an intensity normalization process using a training set of ROIs, followed by a contrast-limited adaptive histogram equalization, and the saturation of the darkest/brightest parts in the ROI were used. After that, a binarization process using Otsu `s thresholding method before keeping the five biggest connected components was applied. After this preliminary step, a DM inspired by Active Shape Models (ASM) [START_REF] Jain | Deformable template models: a review[END_REF] was employed to precisely localize the hippocampus, using a training set of shapes extracted from five to twelve images for each reference slice in the atlas, and deriving a model with eight reference points for the AH and seven points for the DG. This ASM approach is based on an energy minimization framework optimized using a metaheuristic called Differential Evolution (DE) [START_REF] Das | Differential Evolution: A Survey of the State-ofthe-Art[END_REF]. To fully specify the algorithm, we used a crossover rate of 0.9, while the scalar F was set to 0.7. Uniform crossover and DE/target-to-best/1 mutation were employed over a population of 64 individuals and 200 generations. During the registration phase, a classical gradient-based local search method [START_REF] Levenberg | A Method for the Solution of Certain Problems in Least Squares[END_REF] was used in the first step, while Particle Swarm Optimization (PSO) [START_REF] Banks | A review of particle swarm optimization[END_REF] was applied in the second one. Our PSO implementation was run with 24 particles, 40 iterations, c1 = c2 = 2.05, and an inertia factor linearly decreasing with time from 1.0 to 0.1.

The control points of the DM were used to locate different areas of the hippocampus constantly containing both cell bodies and neuropil. They were approximately centered on the region of maximum curvature of the CA1 region and of the CA3 region in the stratum pyramidale of the AH, and in the medial half of the stratum granulosum (sg) for what concerns the DG (Fig. 2A). Afterwards, these areas were first segmented using a system based on Otsu's thresholding method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF], and then further subdivided to obtain 13 different regions of interest (Fig. 2B), able to describe most of the variability of hippocampus' visual features. The 13 regions were then checked to automatically reject incorrect segmentations as follows. A Random Forest (RF) [START_REF] Breiman | Random forests[END_REF] was trained to distinguish between the points lying inside Stratum Pyramidale (sp) and Stratum Granulosum (sg) and those lying outside. Twelve random points were selected in the regions identified in Fig. 1B as CA1 a, CA1 d, CA3 a, CA3 d, DG a, DG d, and classified using the trained RF. Then, the percentage of points correctly classified, along with some statistics of the images containing the regions (area, standard deviation) were processed by another RF, to finally evaluate the segmentation process. If the result was considered unreliable and had to be discarded, the entire process was repeated selecting a different section of the same series from the ABA, and if all slices within the brain region in which the hippocampus is clearly visible had been unsuccessfully processed, the gene was finally discarded. Over 9510 series which were tested, only 255 of were discarded (2.68%), including the ones dismissed due to errors in the ABA, e.g. ruined or missing images. Genes with multiple probes or multiple image series available in the ABA were treated independently because different probes could reveal different transcripts, characterized by dissimilar expression profiles.

Texture Extraction

Textural features of the first and second order were extracted by each region. The size of the windows for each region, as well as the textural features used, are reported in Table 1. The total number of features, i.e. the size of the vector which encodes the visual characteristics of a gene is therefore 220. These features are both of the first order (obtained directly from the image) and of the second order (obtained using Gray Level Co-occurrence Matrix). The Gray Level Co-occurrence Matrix (GLCM) [START_REF] Haralick | Textural Features for Image Classification[END_REF] is a wellestablished method to represent textural information of an image I by defining the distribution of co-occurrence values:

C(i,j)~X m p~1 X m q~1
1 if I(p,q)~i and I(pzDx,qzDy)~j 0 otherwise where m represents the size of the window to analyze, i and j represent the intensity levels of image I (usually the number of levels is reduced to avoid unfeasible sizes of the matrix: in this work The first row contains the 13 regions under consideration (Fig. 2). The second row shows the window sizes used to compute textural features in each region; for example, ''18, 360 means two different sizes (of 18618 and 36636 pixels, respectively). Each feature is the average of the result obtained in windows that cover, as much as possible, the segmented region. doi:10.1371/journal.pone.0074481.t001

we set this number to 16). The offset (Dx, Dy) represents the distance in pixel and the direction between two points. For instance, one can be interested only in vertical or horizontal patterns. In this work, we consider symmetric and non-directional relations between points. Once the matrix is computed, it can be used to extract indices that represent particular features of the image like, for instance, its contrast or its homogeneity.

Model Training and Ranking Procedure

A small dataset composed of 20 genes was used to train the system. The idea was to select a feature subset able both to represent well the dendrite-enriched mRNAs and to successfully distinguish them from negative examples. To do so, we used as positive examples (prototypes) three genes which are well known for the dendrite and spine enrichment of their mRNA, namely Camk2a (ID = 12322) [START_REF] Burgin | In situ hybridization histochemistry of Ca2+/calmodulin-dependent protein kinase in developing rat brain[END_REF][START_REF] Mayford | The 39-untranslated region of CaMKII alpha is a cis-acting signal for the localization and translation of mRNA in dendrites[END_REF], Map2 (ID = 17756) [START_REF] Garner | Selective localization of messenger RNA for cytoskeletal protein MAP2 in dendrites[END_REF] and Arc (ID = 11838) [START_REF] Bramham | Local protein synthesis, actin dynamics, and LTP consolidation[END_REF]. Moreover, we used as negative examples 17 genes with different characteristics. In particular, we chose: Camk2b (ID = 12323), Tubb3 (ID = 22152) and Grin1 (ID = 14810), i.e. three genes whose mRNA is specifically expressed in neurons but is not transported in dendrites; Gapdh (ID = 14433), Pgk1 (ID = 18655) and Pfkm (ID = 18642), i.e. three ubiquitous metabolic enzymes whose mRNA is strongly expressed in neurons but is not transported in dendrites; Gfap (ID = 14580) and Slc1A2 (ID = 20511), i.e. two genes expressed glial cells, particularly in astrocytes; Mag (ID = 17136), Mog (ID = 17441) and Mbp (ID = 17196), i.e. three genes genes expressed in oligodendrocytes; Gad1 (ID = 14415) and Slc6A1 (ID = 232333), i.e. two genes expressed in GABA-ergic interneurons; Slc1A1 (ID = 20510), Slc1A3 (ID = 20512) and Slc17A7 (ID = 72961), i.e. three genes genes expressed in glutamatergic neurons; glial cells, particularly in astrocytes; Sox2ot (320478), a gene producing a non-coding RNA localized in the nucleus. The images chosen for the above genes were reviewed by an expert, who confirmed the correspondence of the expected expression pattern with the pattern revealed by the ABA in situ hybridizations. Features that, in this small dataset, showed a very high correlation between each other (.0.99) were considered as being the same, and one of them was removed. After this, a binary genetic algorithm (with population size 50, run for 300 generations, with crossover rate of 0.8, mutation rate 0.06, and tournament selection of size 4) was used. Every individual in the population encodes a feature subset and the silhouette index [START_REF] Rousseeuw | Silhouettes: a Graphical Aid to the Interpretation and Validation of Cluster Analysis[END_REF] is used to evaluate its effectiveness. This index is computed using the selected subset of features and considering the positive examples in one group and the negative examples in another.

The silhouette index is a measure of an object's similarity to the others of the same group and dissimilarity from the elements in the other groups. In the following formula, ''a'' identifies the average distance (or dissimilarity) between an element with the data in its own group, while ''b'' is the average dissimilarity between the element and the elements of the other group.

s(i)~1 {a(i)=b(i), 0, b(i)=a(i{1), 8 > < > : if a(i)vb(i) if a(i)~b(i) if a(i)wb(i)
The rationale underlying the use of a GA to solve a search space in which, with 220 features, there are 2^220 possible feature subsets, can be found in the ability of such algorithms to solve NPcomplete problems [START_REF] Holland | Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence[END_REF].

Since a genetic algorithm is a stochastic metaheuristic, which may produce different solutions each time it runs, this procedure was repeated 15 times and the features selected in at least 50% of the runs were chosen. This led us to select a subset of 52 features (see Table S1), which were used to process all genes in the ABA and compute the Pearson coefficient between their feature vector and the prototype one, generated by averaging the features of three prototype genes. The results using the 220 features were 0.5066, 0.4302, 0.8150, while using only the selected features were 0.8610, 0.8796, 0.9441; proving the ability of the proposed feature-selection method to facilitate the detection of correlations between genes. It is important to underscore that this pipeline has general applicability and would be able to rank genes according to their similarity to any kind of features, only by providing it a training set including positive and negative examples.

Candidates Validation by qPCR

Total RNA was extracted from brains of adult C57BL/6J mice. All animals were handled in accordance with good animal practice as defined by the relevant national animal welfare bodies, equivalent to the European Convention for the Protection of Vertebrate Animals used for Experimental and other Scientific Purposes (ETS 123). Mouse experimentation protocols were approved by the Nice Sophia Antipolis University regional animal safety committee (CIEPAL-Azur). RNA preparation from total brain and crude synaptosomes was performed as previously described [START_REF] Zongaro | The 3'UTR of FMR1 mRNA is a target of miR-101, miR129-5p and miR-221: implications for the molecular pathology of FXTAS at the synapse[END_REF]. For each sample 500 ng of RNA were retrotranscribed with SuperscriptIII (invitrogen) according to the manufacturer's instruction using both random primers and OligodT. The quantitative real-time PCR was performed on a LightCycler 480 Real-Time PCR System (Roche), using the qPCR Core kit for SYBR Green (Eurogentec), according to the manufacturer's instruction. For each reaction 6.25 ng of cDNA were used. The amplification protocol was: 95uC for 109, followed by 45 cycles of 95uC for 1099, 60uC for 4599 and 72uC for 1099. The relative expression of the transcripts was quantified with the 2 2DCt method [START_REF] Livak | Analysis of relative gene expression data using real-time quantitative PCR and the 2(2Delta Delta C(T)) Method[END_REF] using FMR1 as a reference. For each synaptosomal or brain sample, DCt is the result of the subtraction of the Ct values obtained for Fmr1 (used as a reference, since it is present at equal levels in both synaptosomes and brain) from the Ct values of the non-coding RNAs. This method allows to compare the relative expression levels of a set of transcripts. In all qPCRs a positive control, 2900097C17Rik, was always present. The enrichment of non-coding RNA in the synaptosomal preparation was expressed as the ratio of 2 2DCt between synaptosomal and total brain RNA.

The primers used are: 8030498B09Rik = Forward ATTGGG-TACATGCTCAGGACA, Reverse AGCCAGGGCTACACA-GAGAA; LOC433089 = Forward ATGACCATGGCCTTTT-CATC, Reverse GCTGTGGGGTACAGGGATAA; A830039 N20Rik = Forward CATATCACCCCCGTTGTACC, Reverse TTTTCACTTGGCCAAAAAGC; 2700046G09Rik = Forward CTTGTCCTCTCCTGCACCTC, Reverse AAATAACCAGC-GGGGCTACT; LOC435897 = Forward ATTCCACGTGAT-TGGCAACT, Reverse AAATAACCAGCGGGGCTACT; TC1430156 = Forward TGTCACGGTCAGCTCTGTTC, Reverse AGGGTGGGTCTTCAATTCG; 2900097C17Rik = Forward GACAACGGCCATGTAGTGTG, Reverse ATCC-TATCCCCAAGCCATTT.

Results and Discussion

Detection of Neuropil-enriched Transcripts by the Hippo-ATESC Pipeline

The Hippo-ATESC pipeline was trained using as positive cases three of the best experimentally characterized neuropil-enriched genes described in the literature (Camk2a, Arc and Map2) as well as 17 negative genes (see the methods section for the details). We then set out to scan parasagittal images contained in the ABA to obtain similarity scores between the texture vectors of all informative probes and a prototype, or reference vector, obtained by averaging the vectors of the prototype genes. In particular, to avoid scoring genes expressed at negligible levels in adult mouse hippocampus, we only considered image series in which the expression level or the expression density reported for the hippocampal region or for the hippocampal formation was above 20, a background level determined on the basis of cell cycle genes which are known to be silent in brain after the end of development. For all the genes that resulted above this threshold (n = 9159), we selected parasagittal sections corresponding to levels from 117 to 175 of the reference atlas. In particular, we selected as default the level 145, which we consider as 'center' of the hemisphere, and moved towards the two boundaries if we were not able to get results with the selected slice. For each section we identified the different hippocampal regions and determined the values for the corresponding texture parameters, obtaining vectors of texture features. We then calculated the Pearson correlation coefficients (r) of these vectors with the reference vector. As expected, vectors obtained from sections corresponding to slightly different levels were very similar, as the average Pearson for all the sections of the prototype genes included in the range was equal to 0.8860.17. Of notice, visual inspection of all these sections revealed that the wide range of r-values obtained with the same probe are caused by specific qualitative problems of some sections, which results in very low scores and would therefore produce false negative results, if the procedure should be employed for classification purposes. Nevertheless, genes showing a high Pearson coefficient to the prototype vector displayed a very stable behavior in this test (some examples are shown in Fig. 3). On this basis, we obtained feature vectors for the default level of all ABA probes above threshold and calculated their r values with respect to the reference vector. Interestingly, the distribution of r values did not show inflection points, indicating that it could be naturally used for classification for evidence-based ranking. To evaluate if such an approach could provide valuable new hypothesis for experimental validation, we first analyzed the 20 top scoring probes corresponding to protein-coding genes (Table 2). As expected, a Gene Ontology analysis performed with the DAVID software [START_REF] Dennis | DAVID: Database for Annotation, Visualization, and Integrated Discovery[END_REF] revealed that the most significant common keyword associated to the corresponding genes is ''dendrite'' (P = 0,004). In particular, manual inspection revealed that this list contains some of the best known examples of transcripts localized to dendrites and/or associated with dendritic functions. Indeed, besides to a Camk2a probe not included in the training set, the list contained probes from Dendrin (Dnd) [START_REF] Herb | Prominent dendritic localization in forebrain neurons of a novel mRNA and its product, dendrin[END_REF][START_REF] Kremerskothen | Postsynaptic recruitment of Dendrin depends on both dendritic mRNA transport and synaptic anchoring[END_REF], Psd (also known as Efa6a) [START_REF] Sironi | EFA6A encodes two isoforms with distinct biological activities in neuronal cells[END_REF], microtubule-associated protein 2 (Map2) [START_REF] Garner | Selective localization of messenger RNA for cytoskeletal protein MAP2 in dendrites[END_REF], Git1 [START_REF] Zhang | Synapse formation is regulated by the signaling adaptor GIT1[END_REF] and Spinophilin (Ppp1r9b) [START_REF] Allen | Spinophilin, a novel protein phosphatase 1 binding protein localized to dendritic spines[END_REF][START_REF] Feng | Spinophilin regulates the formation and function of dendritic spines[END_REF]. Accordingly, visual inspection of the corresponding ABA images confirmed a significant signal enrichment in neuropil for many of the probes included in the list Table 2, as it is the case for the Rnf10 gene (Fig. 4). It must be noticed that 7 out the 20 genes are associated to mitochondrial function (Table 2). However, this result is most likely non specific, because it could be easily explained by the high abundance of mitochondria in dendrites. Indeed, the mRNAs of genes mapping to the mitochondrial genome are transcribed within mitochondria, while the mRNAs of nuclear genes encoding mitochondrial oxidative phosphorylation proteins have been found to localize near mitochondria [START_REF] Matsumoto | Localization of mRNAs encoding human mitochondrial oxidative phosphorylation proteins[END_REF]. Therefore, in all these cases a neuropilenriched pattern would be expected. To obtain a more systematic validation, we concentrated on a group of 257 protein-coding gene-derived probes, characterize by a r value of 0.8 or higher.

These probes (Table S2), derived from 117 nuclear genes and from 9 mitochondrial genes, corresponded to 2.7% of all the probes providing in situ hybridization signals significantly higher than background levels. Analysis of the functional annotations of these genes revealed a significant enrichment for phosphoproteins (n = 116, hypergeometric P-value (P) = 2.8E-09) and for GTPbinding proteins (n = 15, P = 4.1E-05), as well as a slight overrepresentation of mitochondrial proteins (n = 30 P = 1.9 E-03). The latter result was even more pronounced if one considers the significant enrichment for probes directly derived from the mitochondrial genome. Besides, it must be considered that other types of ''false positive'' results may be contained in the top-scoring list, such as the Plp1 transcript, which is well known to be expressed by oligodendrocytes [START_REF] Peyron | In situ expression of PLP/DM-20, MBP, and CNP during embryonic and postnatal development of the jimpy mutant and of transgenic mice overexpressing PLP[END_REF].

To evaluate whether the mRNA produced by these genes could be actually enriched in dendrites, we compared our top ranking list with the list of 2250 coding genes that should compose the mouse hippocampal neuropil, according to recent RNA-seq studies [START_REF] Cajigas | The local transcriptome in the synaptic neuropil revealed by deep sequencing and highresolution imaging[END_REF]. This comparison revealed 56 common genes, representing a highly significant intersection (Table S2, P = 2E-07). More importantly, we compared our list with a list of 57 transcript manually annotated by experts as dendritically enriched in the original ABA publication [START_REF] Lein | Genomewide atlas of gene expression in the adult mouse brain[END_REF]. Interestingly, this analysis revealed 12 common genes (P = 2E-17), most of which were concentrated in the top ranks of our list. In comparison, the lists of neuropilenriched transcripts reported in previous micro-dissection-based studies [START_REF] Poon | Identification of process-localized mRNAs from cultured rodent hippocampal neurons[END_REF][START_REF] Zhong | Dendritic mRNAs encode diversified functionalities in hippocampal pyramidal neurons[END_REF] showed a much less significant intersection with the same 'ground truth, because each had in common with it only 2 genes, over a total of 170 and 154, respectively (P-value not significant in both cases). In addition, our list had a significant intersection with those resulting from microdissection studies, having 6 genes in common with the list reported in [START_REF] Zhong | Dendritic mRNAs encode diversified functionalities in hippocampal pyramidal neurons[END_REF] (Pvalue = 0.009) and 9 genes in common with the list reported in [START_REF] Poon | Identification of process-localized mRNAs from cultured rodent hippocampal neurons[END_REF] (n = 154, P = 5E-05). The above numbers are in the same order of magnitude of the intersection between the two latter studies (6 common genes, P = 0.0082). Taken together, these results indicate not only that the Hippo-ATESC pipeline can efficiently score neuropil-enriched transcripts that would be identified by human experts, but also that it can highlight neuropil-enriched transcript which escape visual expert analysis but can be revealed by quantitative methods.

Experimental Validation of Neuropil-enriched Noncoding Transcripts Characterized by High Hippo-ATESC Score

Considering the good results obtained with protein-coding genes, we would expect that the Hippo-ATESC pipeline should be as effective in highlighting neuropil-enriched non-coding RNAs. Therefore, we analyzed the ranking of probes annotated in the ABA as belonging to non-coding genes. The number of probes providing in situ hybridization signals significantly higher than background levels was 414. However, by mapping these probes to a recent version of the mouse genome annotation (UCSC know genes, version mm9), we realized that only 99 of them are not yet associated with protein-coding genes. Only three probes were characterized by a Pearson correlation coefficient .0.8. However, it was very interesting to notice that one of these probes corresponded to the 2900097C17Rik gene, which has been previously identified as strongly enriched in dendrites by human expert-based inspection of ABA [START_REF] Mercer | Specific expression of long noncoding RNAs in the mouse brain[END_REF]. In consideration of the high specificity of this result and of the small number of non-coding RNAs with very high score, we decided to better analyze the other 6 genes characterized by a Pearson correlation coefficient of at least 0.7 (Table 3). We hypothesized that if the RNA encoded by In the SF and MF columns, the genes which were previously associated to dendritic spine or synaptic function and the genes associated to mitochondrial function are marked, respectively. doi:10.1371/journal.pone.0074481.t002

these genes is actually transported in the mouse hippocampal neuropil, it should be enriched in the synaptodendrosomal compartment, as it turns out to be the case for mRNAs encoding proteins of synaptic relevance, for ncRNA or for microRNAs [START_REF] Zongaro | The 3'UTR of FMR1 mRNA is a target of miR-101, miR129-5p and miR-221: implications for the molecular pathology of FXTAS at the synapse[END_REF][START_REF] Subramanian | Gquadruplex RNA structure as a signal for neurite mRNA targeting[END_REF][START_REF] Muslimov | Spatial codes in dendritic BC1 RNA[END_REF][START_REF] Siegel | microRNAs in neurons: manifold regulatory roles at the synapse[END_REF]. Therefore, we tested this hypothesis by analyzing the relative abundance of the corresponding RNA in crude synaptosomes versus total brain, by means of qRT-PCR. As expected, the RNA of the 2900097C17Rik gene is present at significant level in the synaptosomal compartment, comparable to the one of Fmr1 mRNA, that is located at the synapse [START_REF] Schenck | A highly conserved protein family interacting with the fragile X mental retardation protein (FMRP) and displaying selective interactions with FMRP-related proteins FXR1P and FXR2P[END_REF] and whose synaptosome/total brain ratio is equal to 1 (not shown). Conversely, we did not detect significant expression for the 8030498B09Rik, LOC433089 and LOC435897 genes, and we detected low RNA levels for the 2700046G09Rik gene, although in this case the synaptosomal enrichment was comparable to the positive control (Table 3). The A830039N20Rik and TC1430156 RNAs were detected at levels significantly higher than background (Table 3). However, even the TC1430156 RNA was most likely a false positive results, because it did not show significant synaptosomal enrichment and because it corresponds to the sequence of Rian, an imprinted RNA accumulated in nucleus [START_REF] Hatada | Identification of a new imprinted gene, Rian, on mouse chromosome 12 by fluorescent differential display screening[END_REF]. In contrast, the A830039N20Rik gene represents a truly positive result, because its synaptosomal enrichment was even higher than the positive control (Table 3). Accordingly, visual inspection of the ABA sections revealed for these gene a clear granular positivity in the proximal neuropil of pyramidal cells, which is qualitatively comparable to the positivity displayed by the 2900097C17Rik, although it is quantitatively less intense (Fig. 5).

Conclusions

The Hippo-ATESC pipeline displayed a high exploratory ability in the recognition of neuropil-encoded genes from high resolution ABA images. This method can be seen as a data mining tool that can give helpful information to select target genes whose nature will be studied and confirmed using biological tests. Therefore, the list of candidate neuropil-enriched protein-coding genes that we have provided could represent an important resource for the detection of new genes involved in synaptic plasticity. Indeed, by validating some of the results we have identified a new bona-fide ncRNA enriched in synaptosomes. It will be very interesting to address the functional role of this molecule, by altering its expression levels in neuronal cells.

Supporting Information

Table S1 List of the subset of features used for model training. Notes: ''Energy2'' stands for second order feature Energy. The number in brackets represents the size of the window on which the feature was calculated, when more than one window size was used for the region under consideration. (PDF)

Table S2 List of the probes associated to protein-coding genes characterized by a Pearson coefficient of 0.8 or more. Gene names marked by a star correspond to those mapping to the mitochondrial genome. ABA indicates the Allen Brain Atlas experiment number. In the last three columns the gene is marked if it has been identified in the respective study. (PDF)

Figure 1 .Figure 2 .

 12 Figure 1. Schematic representation of the Hippo-ATESC pipeline. doi:10.1371/journal.pone.0074481.g001

Figure 3 .

 3 Figure 3. Boxplot of the correlation values of the feature vectors obtained from 6 slices of 8 randomly selected genes, as compared to the prototype vector. Note the very stable behavior of genes characterized by a high Pearson coefficient. doi:10.1371/journal.pone.0074481.g003

Figure 4 .

 4 Figure 4. In situ hybridization pattern of the indicated proteincoding genes, obtained from the Allen Brain Atlas. Camk2a was included as a positive control while Ccnb1, which encodes for transcripts expressed only in mitotic cells, was included as a negative control to give an idea of background signal levels. doi:10.1371/journal.pone.0074481.g004

Figure 5 .

 5 Figure 5. In situ hybridization pattern of the indicated non-coding RNA genes, obtained from the Allen Brain Atlas. For positive and negative control, see Fig. 4. doi:10.1371/journal.pone.0074481.g005

Table 1

 1

are indicated. doi:10.1371/journal.pone.0074481.g002

Table 1 .

 1 Different features used for the ranking.

	CA1a CA1b	CA1c	CA1d	CA3a	CA3b	CA3c	CA3d	DGa	DGb	DGc	DGd	DGe
	18,36 24	24	14,28	18,36	24	24	14,28	18,36	24	24	14,28	18,36
	First Order	Mean, Standard Deviation, Coefficient of Variation, Skewness, Kurtosis, Energy, Entropy
	Second Order	Contrast, Correlation, Energy, Homogeneity

Table 2 .

 2 List of the top-scoring 20 protein-coding genes defined by the pipeline.

			Entrez
	Gene Name	ABA	gene ID	Pearson SF	MF
	Ddn	71212512	13199	0.964	X
	9530085C10Rik*	73906215	N/A	0.961		X
	Camk2a	79360274	12322	0.960	X
	0610005I03Rik*	74357749	N/A	0.957		X
	Psd	69352896	73728	0.954	X
	Uqcrb	293241	67530	0.954		X
	Mtap2	69549641	17756	0.954	X
	0610009I12Rik*	74357773	N/A	0.953		X
	0610006F12Rik*	74357761	N/A	0.950		X
	Git1	69672880	216963	0.948	X
	Slc25a23	71280844	66972	0.946		X
	Rnf10	294052	50849	0.945
	Muc10	72737	17830	0.943
	Rpl23	70813131	65019	0.939
	7420498E04Rik*	73927706	N/A	0.938		X
	Senp2	227749	75826	0.936
	Bcap29	72739	12033	0.932
	Mt1	67767450	17748	0.930
	Ppp1r9b	68151446	217124	0.930	X
	Tubb2a	69838608	22151	0.928

Gene names marked by a star correspond to genes mapping to the mitochondrial genome. ABA indicates the Allen Brain Atlas experiment number.

Table 3 .

 3 Relative expression in adult brain and synaptosomes of the top-scoring non-coding RNAs identified by the Hippo-ATESC pipeline, as determined by qRT-PCR.In the last column is reported the ratio between synaptosomal and total brain levels. A.U. = arbitrary units. doi:10.1371/journal.pone.0074481.t003

					Synaptosomal levels	Total brain
	Gene Name	ABA	Entrez gene ID	Pearson	(A.U)	levels (A.U.)	Ratio
	8030498B09Rik	70227944	77547	0.873	ND	ND	ND
	2900097C17Rik	71764607	347740	0.856	16.7260.49	17.2161.6	0.97
	LOC433089	71789951	433089	0.831	ND	ND	ND
	A830039N20Rik	69514374	268723	0.776	5.2360.07	3.8360.3	1.37
	2700046G09Rik	69202980	67188	0.764	0.0660.01	0.0660.01	0.93
	LOC435897	71022615	435897	0.744	ND	ND	ND
	TC1430156 (Rian)	74580805	75745	0.708	8.4660.7	13.7160.82	0.62

PLOS ONE | www.plosone.org

September 2013 | Volume 8 | Issue 9 | e74481

Acknowledgments

All the mouse brain images were downloaded from the website: Allen Mouse Brain Atlas. http://mouse.brain-map.org. Seattle (WA): Allen Institute for Brain Science. 2009.

This work was funded by the Neuroscience Program of the ''Compagnia di S.Paolo'' to MG, SC and FDC. FDC is also funded by the Telethon Foundation and by the Je ´rome Lejeune Foundation. Roberto Ugolotti was funded through a fellowship from the same program. Pablo Mesejo was funded by the European Commission (MIBISOC Marie Curie Initial Training Network, FP7 PEOPLE-ITN-2008, GA n.238819). BB is supported by INSERM, CNRS, ANR (E-RARE 09), Fondation pour la Recherche Me ´dicale. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions

Conceived and designed the experiments: RU PM BB MG SC FDC. Performed the experiments: RU PM SZ. Analyzed the data: RU PM SZ GB FB IM. Contributed reagents/materials/analysis tools: RU PM BB SC. Wrote the paper: RU PM SC FDC.