Sandra Bringay
email: bringay@lirmm.fr

Nicolas Béchet
email: nicolas.bechet@inria.fr

Flavien Bouillot
email: bouillot@lirmm.fr

Pascal Poncelet
email: poncelet@lirmm.fr

Mathieu Roche
email: mroche@lirmm.fr

Maguelonne Teisseire
email: maguelonne.teisseire@cemagref.fr

Towards an On-Line Analysis of Tweets Processing

Tweets exchanged over the Internet represent an important source of information, even if their characteristics make them difficult to analyze (a maximum of 140 characters, etc.). In this paper, we define a data warehouse model to analyze large volumes of tweets by proposing measures relevant in the context of knowledge discovery. The use of data warehouses as a tool for the storage and analysis of textual documents is not new but current measures are not well-suited to the specificities of the manipulated data. We also propose a new way for extracting the context of a concept in a hierarchy. Experiments carried out on real data underline the relevance of our proposal.

Introduction

In recent years, the development of social and collaborative Web 2.0 has given users a more active role in collaborative networks. Blogs to share one's diary, RSS news to track the lastest information on a specific topic, and tweets to publish one's actions, are now extremely widespread. Easy to create and manage, these tools are used by Internet users, businesses or other organizations to distribute information about themselves. This data creates unexpected applications in terms of decision-making. Indeed, decision makers can use these large volumes of information as new resources to automatically extract useful information.

Since its introduction in 2006, the Twitter website 5 has developed to such an extent that it is currently ranked as the 10 th most visited site in the world 6 . Twitter is a platform of microblogging. This means that it is a system for sharing information where users can either follow other users who post short messages or be followed themselves. In January 2010, the number of exchanged tweets reached 1.2 billion and more than 40 million tweets are exchanged per day 7 . When a user follows a person, the user receives all messages from this person, and in turn, when that user tweets, all his followers will receive the messages. Tweets are associated with meta-information that cannot be included in messages (e.g., date, location, etc.) or included in the message in the form of tags having a special meaning. Tweets can be represented in a multidimensional way by taking into account all this meta-information as well as associated temporal relations. In this paper, we focus on the datawarehouse [START_REF] Codd | Providing OLAP (On-Line Analytical Processing) to User-Analysts: An IT Mandate[END_REF] as a tool for the storage and analysis of multidimensional and historized data. It thus becomes possible to manipulate a set of indicators (measures) according to different dimensions which may be provided with one or more hierarchies. Associated operators (e.g., Roll-up, Drilldown, etc.) allow an intuitive navigation on different levels of the hierarchy.

This paper deals with different operators to identify trends, the top-k most significant words over a period of time, the most representative of a city or country, for a certain month, in a year, etc. as well as the impact of hierarchies on these operators. We propose an adapted measure, called T F -IDF adaptive , which identifies the most significant words according to level hierarchies of the cube (e.g., on the location dimension). The use of hierarchies to manage words in the tweets enables us to offer a contextualization in order to better understand the content. We illustrate our proposal by using the MeSH8 (Medical Subject Headings) which is used for indexing PubMed articles 9 .

The rest of the paper is organized as follows. Section 2 describes a data model for cubes of tweets and details the proposed measure. In Section 3, we consider that a hierarchy on the words in tweets is available and propose a new approach to contextualize the words in this hierarchy. We present some results of conducted experiments in Section 4. Before concluding by presenting future work, we propose a state-of-the-art in Section 5.

2 What is the most appropriate measure for tweets?

Preliminary Definitions

In this section we introduce the model adapted to a cube of tweets. According to [START_REF] Pérez-Martínez | Contextualizing data warehouses with documents[END_REF], a fact table F is defined on the schema D = {T , . . . , T n , M } where T i (i = 1, .., n) are the dimensions and M stands for a measure. Each dimension T i is defined over a domain D partitioned into a set of categories C j . We thus have D = ∪ j C j . D is also provided with a partial order D to compare the values of the domain D. Each category represents the values associated with a level of granularity. We note e ∈ D to specify that e is a value of the dimension D if there is a category C j ⊆ D such that e ∈ ∪ j C j . Note that two special categories are distinguished and are present on all dimensions: ⊥ D et D ∈ C D corresponding respectively to the level of finer and higher granularity. In our approach, the partial order defined on the domains of the dimensions stands for the inclusion of keywords associated to the values of the dimensions. Thus, let e 1 , e 2 ∈ ∪ j C j be two values, we have e 1 D e 2 if e 1 is logically contained in e 2 .

The data model

We instantiate the data model of the previous section to take into account the different dimensions of description and propose a specific dimension to the words associated to tweets.

Let us consider, for example, the analysis of tweets dedicated to the Duncan diet (e.g., "The Dukan diet is the best diet ever, FACT!!! Its just meat for me for the next 5 day YEESSS"). We wish, for example, to follow the comments or opinons of different people on the effectiveness of a diet. In order to extract the tweets, we query Twitter using a set of seed words: Duncan, diet, slim, protein, etc.. In this case, the original values of the word dimension are dom(word) = { Duncan, diet, slim, protein, ...}.

Figure 1 illustrates the data model. We find the dimension (location ⊥ location = City ≤ State ≤ Country ≤ location ), and the dimension time (⊥ time = day ≤ month ≤ semester ≤ year ≤ time ). The domain of the word dimension is that of the seed words with the words appearing frequently with them. In the fact table, several measures may be used. Traditionally it is the TF-IDF. This issue is addressed in the next section.

Towards an appropriate measure

Relying only on knowledge of the hierarchy in a cube does not always allow a good aggregation (i.e., corresponding to a real situation). For instance, the characteristics of the words in tweets are not necessarily the same in a State and in a City. The aggregation measure that we propose is based on approaches from Information Retrieval (IR).

In our process, the first step is to merge the number of occurrences of words specific to a level. More precisely, we list all the words located in tweets that match a given level (e.g., City, State, Country). If the user wishes to focus the search on a specific City, the words from the tweets coming from this city form a vector. We can apply this same principle to the State by using a Roll-up operator. The aim of our work is to highlight the discriminant words for each level.

Traditionally, T F -IDF measure gives greater weight to the discriminant words of a document [START_REF] Salton | A vector space model for automatic indexing[END_REF]. Like [START_REF] Grabs | ETH Zurich at INEX: Flexible Information Retrieval from XML with PowerDB-XML[END_REF], we propose a measure called T F -IDF adaptative aiming to rank the words according to the level where the user is located and defined as follows:

T Fi,j -IDF k i = ni,j k n k,j × log2 |E k | |{e k j : ti ∈ e k j }| (1)
where |E k | stands for the total number of elements of type k (in our example, k = {City, State, Country}) which corresponds to the level of the hierarchy that the decision maker wants to aggregate. |{ej : ti ∈ ej}| is relative to the number of elements of type k where the term ti appears.

A hierarchy of words for tweets

In this section, we adopt a hierarchy on the words to allow the contextualization of words in tweets.

The data and the model

For the hierarchy, we use the MeSH (Medical Subject Headings) 10 National Library of Medicine's controlled vocabulary thesaurus. It consists of sets of terms naming descriptors in a twelve-level hierarchy that permits the search to be carried out at various levels of specificity. At the most general level of the hierarchical structure are very broad headings such as "Anatomy" or "Mental Disorders". More specific headings are found at more narrow levels, such as "Ankle" and "Conduct Disorder". In 2011, 26,142 descriptors are available in MeSH. There are also over 177,000 entry terms that assist in finding the most appropriate MeSH Heading, for example, "Vitamin C" is an entry term to "Ascorbic Acid". The data model is updated to take into account this new dimension. Compared to the previous model (See Figure 1) the dimension "Word" has been replaced by MeSHWord. MeSHWord has a partial order, M eSHW ord , to compare the different values of the domain. One of the main problems with the use of this thesaurus is that different terms may occur at various levels in the hierarchy. This ambiguity raises the problem of using operators like Roll-up or Drill-down to navigate in the cube. In order to illustrate this problem let us consider the following example.

Example 1 Let us consider the following tweet: "pneumonia & serious nerve problems. can't stand up. possible myasthenia gravis treatable with meds.". If we look in MeSH for the descriptor pneumonia, we find this term occurring in several places (See Figure 2). Depending on the position in the hierarchy, a Rollup operation on pneunomia will not give the same result (i.e., "respiratory tract diseases" versus. "lung diseases").

How to identify the context of a tweet?

We have shown in Example 1, the difficulty of locating the context of a term in the hierarchy. However, a closer look at the tweet shows that the words themselves can be helpful to determine the context. In order to disambiguate polysemous words in the hierarchy of MeSH, we adapt the AcroDef M I 3 method described in [START_REF] Roche | Managing the acronym/expansion identification process for text-mining applications[END_REF] where the authors show the efficiency of this method in a biomedical context. This measure is based on the Cubic Mutual Information [START_REF] Daille | Approche mixte pour l'extraction automatique de terminologie : statistiques lexicales et filtres linguistiques[END_REF] that enhances the impact of frequent co-occurrences of two words in a given context. For a context C, AcroDef M I 3 is defined as follows:

AcroDef C M I 3 (m1, m2) =
(nb(m1 and m2 and C)) 3 nb(m1 and C) × nb(m2 and C)

In our case, we want to calculate the dependence between a word m to disambiguate and different words m t of the tweets using the context of the hierarchy (i.e., parents p of the word m).

Example 2 Let us consider the word 'pneunomia' to disambiguate in Example 1. Here we calculate the dependence between this word m and the other words following 'pneumonia' (nouns, verbs, and adjectives are selected with a Part-of-Speech process): 'serious' and 'nerve'. This dependence is calculated regarding the context of both possible fathers in the MeSH hierarchy. In order to predict where in the MeSH thesaurus we have to associate the word 'pneunomia', we perform the following operations:

nb(pneunomia, mt, "lung diseases") = 227 (number of returned pages with the queries 'pneunomia serious "lung diseases"' and 'pneunomia nerve "lung diseases"') nb(pneunomia, mt, "respiratory tract inf ections") = 496

The dependence of the terms is given by: -AcroDef "lung diseases" M I 3 (pneunomia, mt) = 0.02 -AcroDef "respiratory tract inf ections" M I 3 (pneunomia, mt) = 0.11 Thus, in the tweet from Example 1, for the word pneunomia, we will preferably do the aggregation at the level of the concept 'respiratory tract infections' of the MeSH. Note that this step of disambiguation, which is essential for data from MeSH, is quite costly in terms of the number of queries. It therefore seems more appropriate to call these functions during the ETL process rather than carrying out such processing when browsing the cube.

Experiments

In order to evaluate our approach, several experiments were conducted. These were performed using PostgreSQL 8.4 with the Pentaho Mondrian 3.20 environment. To extract the tweets related to the vocabulary used in MeSH, we focus on the tweets related to "Disease" and queries Twitter by using all the terms of the corresponding hierarchy. We collected 1,801,310 tweets in English from January 2011 to February 2011. In these experiments, we analyze the first words returned by the TF-IDF adaptive (highest scores). For example, the following table presents the first 12 words of tweets in the United States, for the State of Illinois and the City of Chicago during the month of January 2011. Now we consider an example of the application of our approach. Figures 3 and4 visualize the worldwide coverage of the words hepatitis and pneunomia excluding the United States, the United Kingdom, and Canada. This coverage is obtained by fixing the location dimension and by examining the frequency of the Word over the period. Finally we evaluated the prediction measure (i.e., AcroDef M I 3 ) within the MeSH hierarchy (see section 3.2). We extracted more than 5,000 Facebook messages (the same kind of messages as tweets) from the food topic. These messages contain at least one polysemous term (i.e. a term which can be associated to the hierarchy food and beverages) and one or two other hierachies of MeSH: Eukaryota, lipids, plant structures, and so forth. A correct prediction means that AcroDef M I 3 associates this polysemous term with the food and beverages concept. In the following table, three types of elements are used in order to characterize the hierarchy (context of the AcroDef M I 3 measure): Father (F), grand-father (GF), and father + grand-father (FGF). This table shows that (1) the use of more generic information (grand-father) is more relevant, (2) the association of all the available information improves the quality of the prediction. In our future work we plan to add other hierarchical information (e.g. son, cousins).

Elements of the hierarchy used F GF FGF Prediction 60.8% 63.6% 68.0%

Related work

The analysis of textual data from tweets has recently been addressed in many research studies and many proposals exist. For example, in [START_REF] Sakaki | Earthquake Shakes Twitter Users: Real-time Event Detection by Social Sensors[END_REF], the authors propose analyzing the content of the tweets in real time to detect alarms during an earthquake. The authors of TwitterMonitor [START_REF] Mathioudakis | Twittermonitor: trend detection over the twitter stream[END_REF] present a system to automatically extract trends in the stream of tweets. A quite similar approach is proposed in [START_REF] Benhardus | Streaming trend detection in twitter[END_REF]. However, to the best of our knowledge, most existing studies mainly focus on a specific analysis of tweets and do not provide general tools for the decision maker (i.e., for manipulating the information embedded in tweets according to their needs). Thus, few studies have been interested in the use of cubes to manage tweets. Recent work has focused on integrating textual data in data warehouses. In this context, aggregation methods suitable for textual data have been proposed. For example, in [START_REF] Keith | Analyzing large collections of electronic text using olap[END_REF], the authors propose using Natural Language Processing techniques to aggregate the words with the same root or the same lemmas. They also use existing semantic tools such as Wordnet or Roget to group words together. Apart from using morpho-syntactic and semantic knowledge, other studies consider numerical approaches from the field of Information Retrieval (IR) to aggregate textual data. Thus, the authors of [START_REF] Lin | Text Cube: Computing IR Measures for Multidimensional Text Database Analysis[END_REF] propose aggregating documents according to keywords by using a semantic hierarchy of words found in the datawarehouses and some measures from IR. Such methods from IR are also used in the work of [START_REF] Pérez-Martínez | Contextualizing data warehouses with documents[END_REF] which takes into account a "context" and "relevance" dimension to build a datawarehouse of textual data called R-cube. Other approaches add a new specific dimension. For example, in [START_REF] Zhang | Topic cube: Topic modeling for olap on multidimensional text databases[END_REF], the authors add a "topic" dimension and apply the PLSA approach [START_REF] Hofmann | Probabilistic latent semantic analysis[END_REF] to extract the themes representing the documents in this new dimension. Finally, in [START_REF] Pujolle | Fonctions d'agrégation pour l'analyse en ligne (OLAP) de données textuelles[END_REF] the authors propose aggregating parts of documents to provide the decision maker with words specific to the aggregation. In this context, the authors use a first function to select the most significant words using the classical T F -IDF measure.

In this paper we proposed a new approach to analyze tweets from their multidimensional characteristics. The originality of our proposal is to define and manipulate cubes of tweets. We have shown through two different models and applications: no predefined hierarchy on tweets (i.e., diet analysis) and existing hierarchy (i.e., using the MeSH thesaurus), that the analysis of tweets requires the definition of new measures and that a contextualization step is relevant. Future work involves several issues. First we want to extend the proposed approach to take into account opinions or feelings expressed in the tweets. Recent studies analyze the mood of people (e.g., http://twittermood.org/). We want to enhance these approaches by analyzing the content of tweets and thus be able to automatically extract knowledge such as: who are the people who followed a diet and are dissatisfied? Secondly, we wish to consider tweets as available in the form of a stream and propose new techniques for efficiently storing the data.

Fig. 1 .

 1 Fig. 1. The schema related to a diet application

Fig. 2 .

 2 Fig. 2. An example of the MeSH thesaurus

  Colloque : 22nd International Conference on Database and Expert Systems Applications, DEXA 2011, Toulouse, 29 August -2 September 2011 Author-produced version of the article published in Lecture Notes in Computer Science, 2011, 6861(2), 154-161. The original publication is available at http://www.springer.com DOI: 10.1007/978-3-642-23091-2_15

Fig. 3 .Fig. 4 .

 34 Fig. 3. Distribution of the use of the word hepatitis Fig. 4. Distribution of the use of the word pneunomia

http://www.ncbi.nlm.nih.gov/mesh

http://www.ncbi.nlm.nih.gov/PubMed/

Colloque : 22nd International Conference on Database and Expert Systems Applications, DEXA 2011, Toulouse, 29 August -2 September 2011 Author-produced version of the article published in Lecture Notes in Computer Science, 2011, 6861(2), 154-161. The original publication is available at http://www.springer.com DOI: 10.1007/978-3-642-23091-2_15


C Rodriguez-Carranza

P Mukherjee

D Vigneron

J Barkovich

C Studholme

Measurment of folding in surfaces of arbitrary size in human brain development

This paper describes a novel approach to in-vivo measurement of brain surface folding in clinically acquired neonatal MR image data. Specifically it addresses the problem of comparing folding of surfaces of arbitrary area. Most of the current measures of folding are not independent of the area of the surface from which they are derived. Therefore, applying them to whole brains or subregions of different sizes result in differences which may or may not reflect true differences in folding. In this paper we describe alternative approaches to deriving area independent measures. The measures were applied to twelve premature infants (age 28-37 weeks) from which cortical gray and white surfaces were extracted. Experimental results show that previous folding measures are sensitive to the size of the surface patch, and that the area independent measures proposed here provide significant improvements. Such a system provides a tool to allow the study of structural development in the neonatal brain within specific functional subregions, which may be critical in identifying later neurological impairment.

Introduction

The percent of infants born preterm, or at less than 37 completed weeks of gestation, has increased significantly in the last two decades. The latest report [START_REF] Hamilton | Births: Preliminary data for 2004[END_REF] in the USA indicates that almost one in every eight births is preterm. There is growing evidence that premature birth can result in structural and functional alterations of the brain, which are related to adverse neurodevelopmental outcome later in life [START_REF] Hüppi | Structural and neurobehavioral delay in postnatal brain development of preterm infants[END_REF][START_REF] Inder | Abnormal cerebral structure is present at term in premature infants[END_REF]. Some of the challenges that preterm infants face range from spastic motor deficits (cerebral palsy) [START_REF] Volpe | Neurologic outcome of prematurity[END_REF], impaired academic achievement [START_REF] Cooke | Ophtalmic impairment at 7 years of age in children born very preterm[END_REF][START_REF] Hack | Outcomes in young adulthood for very-lowbirth weight infants[END_REF], and behavioral disorders [START_REF] Inder | Defining the nature of the cerebral abnormalities in the premature infant: a qualitative magnetic resonance imaging study[END_REF][START_REF] Luciana | Cognitive development in children born preterm: implications for theories of brain plasticity following early injury[END_REF]. However, the conditions that cause the cerebral abnormalities that underlie these common and serious developmental disabilities are not entirely understood [START_REF] Inder | Abnormal cerebral structure is present at term in premature infants[END_REF]. The wider availability of clinical in vivo magnetic resonance imaging of neonatal brain anatomy, provided by systems which make use of an MRI compatible incubator, creates a new opportunity to quantify brain development. In this work we are particularly interested in the study of the cortical folding, or gyrification, in preterm infants, because it may reflect underlying functional organization.

In our earlier work [START_REF] Rodriguez-Carranza | An iso-surface folding analysis method applied to premature neonatal brain development[END_REF] we applied previously proposed global average measures of folding to surfaces extracted from neonatal MRI, and showed their use in tracking global age changes. The main limitation was that the previously proposed measures were heavily dependent on the size of the surface being examined. Since brain surface area increases dramatically with brain development, previous measures cannot probe whether cortical folding is following a normal pattern, independent of size. Additionally these cannot be calculated on functional sub-regions of the cortex which may have a different surface area in different individuals. To address the problem, in this paper we describe and analyze measures of global folding that are independent of the size of the surface of analysis. We apply the measures to surfaces extracted from MR images of twelve premature infants with ages from 28-37 weeks. The folding of whole surfaces, as well as their left and right hemispheres, was studied.

Theory

Several global 3D brain surface measures of gyrification have previously emerged. Some are based on surface principal curvatures (k 1 , k 2 ), mean curvature (H =

1 2 (k 1 + k 2 )), or Gaussian curvature (K = k 1 k 2 ).
Examples are folding index (F I) [START_REF] Van Essen | Structural and functional analyses of human cerebral cortex using a surface-based atlas[END_REF], the intrinsic curvature index (ICI) [START_REF] Van Essen | Structural and functional analyses of human cerebral cortex using a surface-based atlas[END_REF], L 2 norm of the mean curvature (M LN ) [START_REF] Batchelor | Measures of folding applied to the development of the human fetal brain[END_REF], L 2 norm of Gaussian curvature (GLN ) [START_REF] Batchelor | Measures of folding applied to the development of the human fetal brain[END_REF], and average curvature [START_REF] Magnotta | Quantitative in vivo measurement of gyrification in the human brain: changes associated with aging[END_REF]; global shape index (GS) and global curvedness (GC) were defined in our earlier work [START_REF] Rodriguez-Carranza | An iso-surface folding analysis method applied to premature neonatal brain development[END_REF] based on the local shape descriptors curvedness (c) [START_REF] Koenderink | Surface shape and curvature scales[END_REF] and shape index (s) [START_REF] Koenderink | Surface shape and curvature scales[END_REF]. Examples of non-curvature based measures are the gyrification index [START_REF] Zilles | The human pattern of gyrification in the cerebral cortex[END_REF], calculated as the ratio of the entire cortical contour of the brain to the superficially exposed contour, and roundness1 (Rn), based on surface area (A) and volume (V ). A list of current global 3D measures of folding is shown on the left column of Table 1. These expressions are normalized to yield the unit value for a sphere [START_REF] Rodriguez-Carranza | An iso-surface folding analysis method applied to premature neonatal brain development[END_REF].

Size-independent measures of surface folding

Examination of the form of the current global measures in Table 1 reveals a critical dependency on surface area, with the exception of global shape index and gyrification index. This can be better illustrated with the following example, using M LN to measure folding. Take a whole sphere of radius Ro, half a sphere of radius Ro, and a whole sphere of half that radius. The three objects have the same surface complexity, hence a measure of global folding should yield identical results. By definition, H = 1/R at each point in a sphere. The whole sphere with radius Ro yields

M LN = 1 4π A H 2 = 1 4π × 4πRo 2 × (1/Ro) 2 = 1
, and so does the whole sphere with radius Ro/2: M LN = 1 4π × (4π Ro 2 ) 2 × (2/Ro) 2 = 1. On the other hand, M LN for the half sphere of radius Ro is half that value:

M LN = 1 4π × 1 2 4πRo 2 × (1/Ro) 2 = 1 2 .
This demonstrates the dependency of the measure on the size of the surface of analysis. Normalizing M LN with the surface area does not alleviate the problem:

M LN = 1 A A H 2 = 1 4πRo 2 (4πRo 2 )(1/Ro) 2 = 1
Ro 2 (whole sphere of radius Ro). We propose two normalization factors: 1) T = 3V /A and 2) |H| = 1

A | A H|. The resulting T -normalized and H-normalized expressions (T -measures and H-measures, respectively) are shown in the right column of Table 1. For a sphere, T = Ro, |H| = 1/Ro, and it can be verified that M LN T = M LN H = 1 for the three sphere cases. Similarly for the rest of the newly defined area-independent measures.

M LN = 1 4π A H 2 [1] M LN T = T 2 A A H 2 GLN = 1 4π A * A K 2 [1] GLN T =T 4 1 2 A K 2 ICI= 1 4π A K + [4] ICI T =T 1 A A K + F I= 1 4π A ak [4] F IT =T 1 A A ak GC= 1 √ A4π A c [16] GC T =T 1 A A c U = 1 A A U [13] U T = T A U A U V = 1 A A V VT =T 1 A V A V Rn= A 3 √ 36πV 2 [1] SH2SH=T A H 2 / A H SK2SK=T A K 2 / A K Area-independent GS= 1 A A s [16] M LNH = 1 |H| A H 2 /A Gyrification Index [21] GC H = 1 A|H| A c AF i =A i /A, i ∈ {H + , K + } Notation: c = 1 2 (k 2 1 + k 2 2 ) (curvedness [10]), s = 2 π arctan k 2 +k 1 k 2 -k 1 (shape index [10]), ak = |k 1 |(|k 1 | -|k 2 |), T = 3 V A , and |H| = 1 A | A H|. U ∈ [H + , H -], V ∈ [K + , K -],
and A H+ (or A K+ ) is the area of the surface with positive H (or K) curvature.

In addition, we defined three new area-independent folding measures: SH2SH, SK2SK, and AF {H + ,K -} . The first two were based on the rationale that measures composed of ratios of local curvature factors would intrinsically eliminate the dependence on area. The latter was based on the idea that a reasonable characterization of the degree of folding in a brain surface is the fraction of the surface which contains convex folds (gyri). This can be mathematically characterized by looking at the relative portion of the surface which has positive mean or Gaussian curvature. For an undeveloped brain, the fraction would be large, and it would progressively decrease as the concave folds (sulci) appear. The definition of the three expressions is shown in Table 1.

Method

The sphere example was useful to identify global folding measures that have the potential to consistently evaluate surface folding, independent of the size of surface of analysis. A more formal verification of this area-independence property is in progress. In this work we present an experimental comparison of the measures on two datasets. The first dataset consisted of 15 neonatal brains. The surfaces of this brains not only increase in size with age, but they also become more complex. So, folding varies proportionally to size. The second dataset consisted of iso-surfaces for which surface folding varied inversely to size. They were generated from one neonatal brain at various percentage occupancy thresholds (see Section 3.3). A larger threshold produces a smaller brain and larger separation between sulci walls, which translates to increased folding. For both datasets, gyrification was measured on left and right hemispheres, in addition to the whole brains. The infants in this study are normal, therefore gyrification is expected to be similar on both sides. Area-independent measures are expected to yield similar values for each hemisphere and the full-brain, but not so the area-dependent measures.

Data

Human Neonatal Images. High resolution (0.703 × 0.703 × [1.5 -2.2] mm) 3D T1 weighted SPGR images were acquired on premature infants using a 1.5T GE MRI scanner with an MRI compatible incubator. The gestational ages (GA) of the twelve premature infants in this study ranged between 24-31 weeks. The postmenstrual ages (gestational age + postnatal age) at the time of acquisition were 28-37 weeks. A subsequent scan for three infants was available, hence a total of 15 brains were processed.

Image Segmentation

The outer gray matter and the gray-white matter interface2 surfaces of the premature brains were extracted semi-automatically as described in [START_REF] Rodriguez-Carranza | An iso-surface folding analysis method applied to premature neonatal brain development[END_REF] using the rview software package [START_REF] Studholme | [END_REF]. Cortical gray matter was segmented for fifteen brains and gray-white matter interface for thirteen brains. A value of 1000 was assigned to brain voxels and 0 to background voxels.

Computation of Brain Surface Folding

Each binary segmented brain was supersampled using voxel replication to prevent loss of fine scale features (2 × 2 × 4). A voxelwize approach to surface curvature estimation from iso-surfaces, derived from that developed by Rieger [START_REF] Rieger | On curvature estimation of iso surfaces in 3D gray-value images and the computation of shape descriptors[END_REF], was then employed. It avoids the need of a parametric model. A summary of the sequence of steps (described in detail in [START_REF] Rodriguez-Carranza | An iso-surface folding analysis method applied to premature neonatal brain development[END_REF]) is as follows: 1) computation of the image gradient g = ∇f (x, y, z) of the replicated binary volume; 2) computation of the gradient structure tensor (GST) [START_REF] Rieger | On curvature estimation of iso surfaces in 3D gray-value images and the computation of shape descriptors[END_REF] defined as T = gg t ; 3) calculation of the eigenvectors v 1 , v 2 , v 3 of T and the mapping M (v 1 ) =

v 1 v t 1 ||v1|| ; 4) on the 50% occupancy iso-surface, computation of the principal curvatures from |k 1,2 | = 1 √ 2 ||∇ v 2,3 M (v 1 )
|| F (Fröbenius norm); the sign of k 1 and k 2 is determined from the Hessian matrix of the image and the eigenvectors v 2 and v 3 . The curvatures k 1 and k 2 were then used to compute mean (H) and Gaussian (K) curvature and the folding measures of Table 1. Each step involves image smoothing or differentiation. For all images, the 3D filters used at the ith step had a full-width at half maximum (f i ) of:

f 1 = 2.1mm, f 2 = 1mm, f 3 = 2mm, f 4 = 1mm.
The rationale behind the choice of values was to create a smooth surface that would yield smooth variations in curvature. Normalized Gaussian derivatives [START_REF] Lindeberg | On scale selection for differential operators[END_REF] were used for all computations.

The input to this algorithm is the binary replicated segmented brain, for which foreground voxels have a value of 1000 and background voxels a value of 0. Near the border of the brain the kernel incorporates brain and background voxels, hence the raw effect of the convolution on the image is that the resulting volume has voxel values in the range [0,1000]. The resulting voxel values are interpreted as the partial occupancy by brain tissue. The larger the number indicates that more foreground voxels were present in the kernel. For each brain, the surface on which curvature and the global measures of folding were calculated was taken to be the iso-surface with at least 50% occupancy (i.e. voxel values of at least 500) and satisfying 6-connectivity with the background.

Cortical Partitioning using Spatial normalization

Full-brains and their left and right hemispheres were processed with the algorithm just described. The procedure used to identify the left and right hemispheres in each neonatal brain was the following. A reference neonatal anatomy was manually partitioned into left and right hemisphere. A multi-resolution B-Spline based spatial normalization, adapted from that described in [START_REF] Studholme | Accurate template-based correction of brain MRI intensity distortion with application to dementia and aging[END_REF][START_REF] Studholme | Deformation tensor morphometry of semantic dementia with quantitative validation[END_REF] and as used in [START_REF] Cardenas | Using automated morphometry to detect association between ERP latency and structural brain MRI in normal adults[END_REF], was then applied to estimate a spatial transformation mapping from this reference to each subject MRI being studied. This transformation was then numerically inverted to allow the assignment of the nearest reference label to a given subject image voxel, which brought the voxelwize partitioning of the reference anatomy to the same space as the surface extracted from each infant. The partitioning was then used to constrain the voxelwize evaluation of the folding measures.

We applied the measures to surfaces from the 15 MRI studies to investigate their response to folding with age. The relationship between age and folding was analyzed with regression and for this work a linear model was used. Due to space limitations results from only seven measures are shown in Figure 1. Top two rows correspond to gray matter results, and bottom two rows to white matter. The first observation is that measures with large goodness-of-fit for gray matter, had very small goodness-of-fit for white matter, and viceversa. This was true for the whole brain and both hemispheres. The same was observed in our results [START_REF] Rodriguez-Carranza | A system for measuring regional surface folding of the neonatal brain from MRI[END_REF] on 10 MRI studies. The measures with greater goodness-of-fit for gray matter were the T -normalized measures: ICI T , F I T , H + , K -, K + ; for white matter, GS (global shape index) and the |H|-normalized measures: M LN H and GC H . The goodness-of-fit for gray matter and white matter (in that order) on whole surfaces (depicted as * in the figure) was: CI T (0.85, 0.004), F I T (0.85, 0.19), GLN T (0.80, 0.21), ICI T (0.83, 0.002), M LN T (0.80, 0.07), GS (0.61, 0.95), SK2SK (0.67, 0.46), SH2SH (0.77, 0.25), M LN H (0.48, 0.96), GC H (0.52, 0.95), H -(0.73, 0.07), H + (0.85, 0.12), K -(0.85, 0.27), K + (same as ICI T ), AF H + =(0.62, 0.87), AF K + =(0.45, 0.13). The measures had similar goodness-offit scores for the left and right hemispheres. Area, volume, and T varied linearly in all cases.

Next we analyzed the dataset for which smaller brains were more folded. The iso-surfaces were created from a single brain (age 28 weeks) at percentage occupancy thresholds in the range 5-95% (in increments of 5%). An example of one coronal slice from five iso-surfaces obtained at five different thresholds is shown in Figure 2. The results for five measures on whole surfaces and their left and right hemispheres are shown in Figure 3. In contrast with the first dataset, the folding measures had similar goodness-of-fit for both gray and white matter (either both large or both small). Another difference was that F I T presented a small goodness-of-fit. For this iso-surface dataset the area of gray matter varied in a non-linear fashion (plot not shown). This behavior is explained by the separation of the gray matter volume into two hemispheres, for larger thresholds (see Figure 2). The increase in surface area due to this separation is dominant over the relative decrease in area due to the larger threshold.

It can be observed that in the first dataset (Figure 1) the slopes of Tnormalized measures have different tendencies for gray matter (positive) and white matter (negative). This difference is explained by the slope of T itself. For gray matter the rate of change of volume, with age, was faster than that of area, therefore T = 3V /A had a positive slope. The inverse occurs for white matter. For the second dataset (Figure 3), the slope of T was negative for both gray and white matter, therefore the slope of all T -normalized measures is also negative. Once this slope difference is considered, it is clear that the new measures (Tand |H|-normalized) change consistently with folding, independent of the size of the surface.

Results on both datasets for previously defined measures corroborate their dependency on the size of the surface of analysis. For each hemisphere, the Plots of seven folding measures applied to gray matter (top two rows) and white matter (bottom two rows). Filled shapes correspond to results for left hemisphere, and empty ones for right hemisphere. Results for the whole surface ( ) are also included. magnitude of previously defined measures was at least half the value than for the whote-brain. In contrast, for the newly defined measures, their magnitude for whole-brain was in between that for left hemisphere and right hemisphere. As an example, results for M LN and GLN are shown in Figures 1 and3, respectively.

Assessment of intra-observer segmentation variability (preliminary stage). The white matter of two brains, from the 15 available, was segmented and analyzed with an approach similar to the one in [START_REF] Batchelor | Measures of folding applied to the development of the human fetal brain[END_REF]. One brain was 29.1 weeks and was segmented four times; the other brain was 32.4 weeks and was segmented twice. Measures were computed on all the segmented surfaces. Given a set of segmentations for one brain, the segmentation variability was assessed as the maximum percentage variation of a measure from the mean value in the set. This was compared to the percentage change of the mean value of a measure between the 29.1 week brain and the 32.4 week brain. This number was taken as an indication of the extent of the variability in segmentations. An example of the results is given for M LN H , GS, and AF H + . In that order, the maximum percentage variation for each measure on the 29.1 week old brain was: 2.81%, Fig. 2. A coronal slice each from five gray matter (top) and white matter (bottom) iso-surfaces extracted at percentage occupancy thresholds: 10%,30%,50%,70%,90%. All slices come from one brain. 3. Plots of five folding measures applied to iso-surfaces of a single brain at various occupancy thresholds. Top row shows results for gray matter and bottom row for white matter. Filled shapes correspond to left hemisphere and empty ones to right hemisphere. Results for the whole surface ( ) are also included.

2.57%, and 0.55%. The percentage difference between the 29.1 week brain and the 32.4 week brain was: 35.4%, 41.9%, 8.1%.

Discussion

An understanding of the cortical folding process in the development of premature infants may be important in explaining and predicting abnormal neurological outcome. The use of formal mathematical descriptions provides a more quantitative tool to study the folding process than is available with simple visual evaluation of MRI scans. The long term goal of our research is to create a model to track preterm in vivo neonatal brain cortical development that will help characterize normal gyrification and departures from it.

In this work we have shown that most previously proposed folding measures [START_REF] Batchelor | Measures of folding applied to the development of the human fetal brain[END_REF][START_REF] Van Essen | Structural and functional analyses of human cerebral cortex using a surface-based atlas[END_REF] are dependent of the size of the surface patch on which they are calculated. To alleviate this problem we examined two approaches to area normalization of the measures and proposed new forms of measure, which satisfy the requirement of area-independence. These were evaluated on neonatal brain surfaces. We evaluated the measures ability to detect change in normal development by computing the linear regression of each measure with age. Finding the relationship between folding measures and age is an important part of the goal of understanding brain development. We will explore higher order models when more data is available.

In the dataset of 15 neonatal brains, folding varied proportionally to size. White matter folding tracked best with age, while no single measure rated consistently well for both gray and white matter surfaces. For the iso-surfaces dataset, where folding varied inversely to size, most measures rated well, and similarly for both gray and white matter. Despite differences in the behaviour of the normalization factor T and in the linearity of area on both datasets, most newly defined measures appeared to consistently assess folding. Recommended measures for gray matter are ICI T , H + , K -; and for white matter GS, M LN H , GC H , and AF H + . The proposed normalization factors T and |H| opened the possibility to adapt previously defined measures to quantify gyrification in subregions of the brain. The slope of T appears to determine that of the T -normalized measures. Awareness of this is important for the correct interpretation of results. Alternative normalization factors need to be explored.

Preliminary results indicate that the effect of segmentation variability on the measures seems to be minimal: the differences in the measures with segmentation seems to be at least 10 times smaller than the differences with gestational age. We need to further investigate the sensitivity and validity of the measures in a larger cohort, containing both normal and pathological cases. In conclusion, the proposed new and normalized measures provide area-independent assessment of folding which provide the ability to study local gyrification. Such measures will provide a new tool in assessing global and local surface folding which is independent of the overall surface area, and are hence applicable to developing a model that tracks development in premature infants

  Fig.1. Plots of seven folding measures applied to gray matter (top two rows) and white matter (bottom two rows). Filled shapes correspond to results for left hemisphere, and empty ones for right hemisphere. Results for the whole surface ( ) are also included.

  Fig.3. Plots of five folding measures applied to iso-surfaces of a single brain at various occupancy thresholds. Top row shows results for gray matter and bottom row for white matter. Filled shapes correspond to left hemisphere and empty ones to right hemisphere. Results for the whole surface ( ) are also included.

Table 1 .

 1 Measures of surface folding

	Current measures	New measures
	Area-dependent	Area-independent

The terms surface complexity[START_REF] Magnotta | Quantitative in vivo measurement of gyrification in the human brain: changes associated with aging[END_REF] and isoperimetric ratio[START_REF] Batchelor | Measures of folding applied to the development of the human fetal brain[END_REF] have also been used.Mathematical Foundations of Computational Anatomy (MFCA'06)

In the remaining of the text we will refer to the gray-white matter interface simply as white matter.

Mathematical Foundations of Computational Anatomy (MFCA'06)

Acknowledgments

This work was funded by the NIH grant R01 MH65392.


A D Leow

M C Chiang

S C Huang

A W Toga

P M Thompson

Realizing Unbiased Deformation: A Theoretical Consideration

Keywords: Mutual information, Image registration, Computational anatomy

Maps of local tissue compression or expansion are often recovered by comparing MRI scans using nonlinear registration techniques. The resulting changes can be analyzed using tensor-based morphometry (TBM) to make inferences about anatomical differences. Numerous deformation techniques have been developed, although there has not been much theoretical development examining the mathematical/statistical validity of each technique. In this paper, we propose a basic principle that any registration technique should satisfy: realizing unbiased test statistics under null distribution of the displacement. In other words, any registration technique should recover zero change in the test statistic when comparing two images differing only in noise. Based on this principle, we propose a fundamental framework for the construction and analysis of image deformation. Moreover, we argue that logarithmic transform is instrumental in the analysis of deformation maps. Combined with the proposed framework, this leads to a theoretical connection between image registration and other branches of applied mathematics including information theory and grid generation.

Introduction

N on-linear image registration is a well-established field in medical imaging with many applications in functional and anatomic brain mapping, image-guided surgery, and multimodality image fusion [START_REF] Thompson | A framework for computational anatomy[END_REF][START_REF] Grenander | Computational anatomy: An emerging discipline[END_REF][START_REF] Chung | A unified statistical approach to deformation-based morphometry[END_REF][START_REF] Avants | Geodesic estimation for large deformation anatomical shape averaging and interpolation[END_REF][START_REF] Shen | Very high-resolution morphometry using mass-preserving deformations and HAMMER elastic registration[END_REF][START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF][START_REF] Miller | Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms[END_REF][START_REF] Collins | Automated 3D nonlinear deformation procedure for determination of gross morphometric variability in human brain[END_REF]. The goal of image registration is to align, or spatially normalize, one image to another. In multi-subject studies, this serves to reduce subject-specific anatomic differences by deforming individual images onto a population average brain template.

The deformations that map each anatomy onto a common template can be analyzed voxel-wise to make inferences about relative volume differences between the individuals and the template or statistical differences in anatomy between populations [START_REF] Studholme | Detecting spatially consistent structural differences in Alzheimer's and fronto-temporal dementia using deformation morphometry[END_REF][START_REF] Thompson | Mapping Cortical Change in Alzheimer's Disease, Brain Development, and Schizophrenia[END_REF][START_REF] Thompson | Growth patterns in the developing brain detected by using continuum mechanical tensor maps[END_REF][START_REF] Cachier | Symmetrization of the Non-Rigid Registration Problem using Inversion-Invariant Energies: Application to Multiple Sclerosis[END_REF][START_REF] Leow | Asymmetrical Effects of Lithium on Brain Structure Mapped in Healthy Individuals[END_REF]. Similarly, in longitudinal studies it is possible to visualize structural brain changes that occur over time by deforming subjects' baseline scans onto their subsequent scans, using the deformation map to quantify local changes. This area of computational anatomy is known as tensor-based morphometry [START_REF] Thompson | Growth patterns in the developing brain detected by using continuum mechanical tensor maps[END_REF][START_REF] Cachier | Symmetrization of the Non-Rigid Registration Problem using Inversion-Invariant Energies: Application to Multiple Sclerosis[END_REF][START_REF] Leow | Asymmetrical Effects of Lithium on Brain Structure Mapped in Healthy Individuals[END_REF][START_REF] Leow | Longitudinal Stability of MRI for Mapping Brain Change using Tensor-Based Morphometry[END_REF][START_REF] Ashburner | High-dimensional image registration using symmetric priors[END_REF]. As there are numerous techniques for non-linear registration, one may ask whether some approaches are better than others in practice and/or in theory. In the following section, we will try to answer this fundamental question from a statistical standpoint.

Criteria for proper image registration

One could not study non-linear registration without closely examining the common practice of applying logarithmic transformation to Jacobian maps. Log-transformation of a Jacobian determinant field has become standard practice in most tensor-based morphometry (TBM) papers [START_REF] Thompson | Growth patterns in the developing brain detected by using continuum mechanical tensor maps[END_REF][START_REF] Cachier | Symmetrization of the Non-Rigid Registration Problem using Inversion-Invariant Energies: Application to Multiple Sclerosis[END_REF]. The Jacobian determinant of a diffeomorphic (smooth) map is bounded below by zero but unbounded above. Thus, at any voxel, its null distribution would be a better fit to a symmetric Normal distribution if the Jacobians are logged. Additional arguments that favor logtransformation of Jacobian maps come from the need to symmetrize the probabilities of expansions and shrinkages that are reciprocals of each other. We will discuss two different, yet related, concepts of symmetrizing; each supports the use of log transformation.

Two types of symmetry favor the log transformation of Jacobian maps

The first symmetry operates on a voxel level. The motivation behind this approach is the result of voxel-wise statistical tests. Consider, for example, testing for the presence of voxel-wise mean structural change in a multi-subject experiment. One might want to employ classical statistical approaches and test the null hypothesis of zero mean change using one-sample Student's t tests. A more general linear model typically relies on the residuals being a good approximation to a Gaussian random field. This null model assumes that at each voxel the observed rates of change over time (or relative volumes in a cross-sectional study) collected from all subjects follow a normal distribution, thus suggesting the use of log-transformotherwise the Jacobian determinants are bounded below by zero but unbounded above.

The second type of symmetry is the symmetry of Jacobian distributions inside a region (in contrast, the symmetry previously described deals with multiple observations across subjects, at a single voxel). This second, albeit less intuitive, level of symmetry sheds light on how one might construct better non-linear image registrations.

In theory, a proper registration algorithm should produce unbiased estimators of the real anatomical change. An unbiased algorithm should detect no statistically significant change between two serial images if there is no difference other than noise, and statistics quantifying change collected from all voxels should ideally have a zero mean value; rates of change deviating from zero should be considered as errors. In a classical statistical setting, one might expect that rates of change thus follow a Gaussian distribution with zero mean, again justifying the use of the logarithmic transform.

These two types of symmetry operate at different levels: one across all voxels when comparing two images, and the other across all subjects at one single voxel. Assuming this symmetric lognormal distribution for the residuals in a statistical model for the Jacobian determinant, one symmetrizes the rate of change distribution by considering a halving or doubling of volume to be equally likely a priori. A related approach is taken by Pennec et al. [START_REF] Pennec | Riemannian elasticity: a statistical regularization framework for non-linear registration[END_REF], where the Cauchy-Green strain tensor [START_REF] Arsigny | Fast and simple calculus on tensors in the Log-Euclidean framework[END_REF][START_REF] Woods | Characterizing volume and surface deformations in an atlas framework: theory, applications, and implementation[END_REF] of a deformation mapping is logged and used as a term in a penalty functional that is integrated over the image domain to regularize the deformation.

Detecting no change in the absence of real change as a necessary criterion for proper registrations: realizing unbiased test statistics

As noted previously, the logarithmic transform is fundamental in analyzing deformation maps through Jacobian maps, and an ideal registration algorithm should yield Jacobian maps and test statistics that imply zero-change, when no difference is present between two registered images. We refer to this as the principle of realizing unbiased test statistics under null distribution of the displacement.

Before we show how to construct unbiased test statistics, we first define both T and S, on an image domain Ω, as the two images to be registered. Let us also assume, without loss of generality, that the volume of this domain is 1, i.e., | Ω |= 1. We seek to estimate a transformation h such that S is non-linearly registered to T when deformed by h (i.e., S(h(x)) ). In this paper, we will restrict this mapping to be differentiable, one-to-one, and onto from the image domain onto itself [START_REF] Leow | Inverse Consistent Mapping in 3D Deformable Image Registration: Its Construction and Statistical Properties[END_REF] (in practice/implementation, the one-to-one and onto property can be approximated by extending the boundary towards infinity). Let us associate three probability density functions, defined on this image domain, to the identity mapping (id) as well as the deformation h and its inverse:

pdf h (x) = Dh(x) ; pdf h -1 (x) = Dh -1 (x) ; pdf id (x) = 1.
(1)

Here, the Jacobian matrix of a transformation h is denoted Dh. As noted in [START_REF] Leow | Inverse Consistent Mapping in 3D Deformable Image Registration: Its Construction and Statistical Properties[END_REF], integrating the log of the Jacobian determinant over the image domain simply calculates the Kullback-Leibler distance between the identity map and h.

log Dh(x) dx

Ω ∫ = -KL( pdf id , pdf h ) ≤ 0. (2)
Here KL, the non-negative asymmetric Kullback-Leibler (K-L) distance, between two PDFs p and q, is defined as

KL( p,q) = plog p q dx Ω ∫ ≥ 0;
Given Eq. ( 2), one might ask if we could use the negative integral of the logged Jacobian determinantlog Dh(x) dx Ω ∫ as a basis for regularization. This integral only evaluates to zero when h is volume-preserving everywhere. Minimizing this integral gives us a volumepreserving map, and thus realizes unbiased statistics. On a global scale, this integral also evaluates the overall mean log-Jacobian with respect to the computation domain. Moreover, in the field of information theory, the KL distance is geometrically important, providing a means to compare probability density functions on general manifolds [START_REF] Cencov | Statistical Decision Rules and Optimal Inference[END_REF]. Interestingly, the integral Eq.( 2) has skew-symmetry with respect to h and its inverse

KL( pdf id , pdf h ) = - ∫ log Dh dx = KL( pdf h -1 , pdf id ) = Dh -1 log Dh -1 dx ∫ (3)
Notably, the idea of integrating the square root of the Jacobian, as proposed in [START_REF] Nielsen | Statistical warps, a least committed model[END_REF], to remove the skew-symmetry, is equivalent to calculating the Bhattacharyya distance B, a well-known measure in information theory [START_REF] Nielsen | Statistical warps, a least committed model[END_REF][START_REF] Bhattacharyya | On a Measure of Divergence Between Two Statistical Populations Defined by their Probability Distributions[END_REF].

B( pdf id , pdf h ) = Dh(x) 1/ 2 dx ∫ = Dh -1 (x) 1/ 2 dx = ∫ B( pdf h -1 , pdf id ) (4)
Surprisingly, notice that here the Bhattacharyya distance is symmetrical with respect to its two arguments, as well as inverse-consistent. To connect the KL-distance and Bhattacharyya distance, one can consider the geodesic linking of the two PDFs: P(.,t), parameterized by time t:

P(x,t) = pdf id (x) t pdf h (x) 1-t N , N = pdf id t (x) pdf h (x) 1-t dx ∫ (5)
The Bhattacharyya distance corresponds to the arbitrary choice of t=1/2, while a generalization of the above leads to the Chernoff distance in information theory [START_REF] Cover | Elements of Information Theory[END_REF].

Realizing unbiased deformation using symmetrical KL distance

To construct unbiased deformations, we first generalize eq. ( 2): given any diffeomorphism g mapping Region of interests S (with size a) to T (with size b), we have the

following log | Dg(x) | dx S ∫
≤ logb /a, while equality holds if and only if the Jacobian of g takes the constant value b/a. The implication of this generalization is that, assuming ROI mapping from S to T, the Log operation is unbiased if and only if the corresponding Jacobian field is evenly distributed. Thus, minimizing the negative integral of logged Jacobian (treated as a cost function) as in 2.2 again leads to unbiased statistics. Given this generalization, let us now introduce the following by combining Eq.(3), and its counterpart contributed by the inverse mapping:

( ) ( )

dx x Dh x Dh dx x Dh x Dh pdf pdf KL pdf pdf KL pdf pdf KL pdf pdf KL pdf pdf KL pdf pdf KL pdf pdf KL pdf pdf KL id h h id h id h id h id id h id h id h ) ( log 1 ) ( ) ( log 1 ) ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( 1 1
1 1 1 1 - - ∫ ∫ - = - = + = + = + = + - - - - (6)
In this symmetric form, the integrand Dh -1 (

) log Dh is always non-negative, compared to the integrand in Eq. (3) where locally negative numbers can be obtained (when the Jacobian is greater than zero), though globally non-negative. The unbiased property has two layers of meaning: the first one being realizing null distribution when comparing identical images differing in noise (from Eq. ( 2)); the second being realizing correct log Jacobian statistics in ROI analysis (from the generalization of Eq. ( 2)). Under this framework, constructing deformations can be viewed as quantifying the symmetric KL distance between the identity map and the resulting deformation (or the inverse deformation due to its symmetry). Moreover, this framework embeds statistical analyses into the construction of deformations, penalizing deformations that skew the distribution of test statistics. A second interpretation of Eq.( 6) is that it simply calculates the mean log Jacobian for h and its inverse inside the domain, thus encoding regional volume changes. Let us also point out an interesting observation by applying the square root Jacobian integrand as in Eq.( 4) to the log-transformed Jacobian framework. A simple change of variables verifies the following: log Dh

( ) ∫ Dh 1/ 2 dx + log Dh -1 ( ) ∫ Dh -1 1/ 2 dx = 0 (7)
To further link Eq.( 6) to other branches of mathematics, optimization problems involving Jacobian operator are commonly encountered in grid generation [START_REF] Liseikin | Grid Generation Methods[END_REF] and in continuum mechanics, where the Hencky tensor arises as logged tensor parameters in modeling very large deformations. However, we believe that the logarithmic transform has not been formally introduced in the grid generation literature and may also be useful there.

Logarithmic transform: a ubiquitous operation

With Eq. ( 6), realizing unbiased log-Jacobian values can be thought of as equivalent to minimizing the KL-distance. Moreover, in this case the concept of inverse-consistency translates to the symmetrization of the KL-divergence. The choice of the logarithmic operation in defining KL-divergence, though arbitrary at first sight, can now be easily justified by linking it to the statistical analyses of deformations in non-linear registration.

Finally, we comment briefly on how this framework can be implemented. Given an image matching cost function C, we seek, among all deformations minimizing this matching cost, the deformation with minimal symmetric KL distance as defined in Eq [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF]. In practice, this often means implementing a Lagrange multiplier, resulting in the following combined minimization problem

1 h id id h H h pdf pdf KL pdf pdf KL h S S h T T C + + - ∈ λ o o (8) ( ) ) , ( ) , ( ) , , , ( min arg
Here λ is the multiplier, and H is the solution space, of which a common choice consists of all one-to-one, onto, and differentiable maps. Often, the solution is numerically obtained by recursive smoothing or regularization applied to the force field. In this paper, we focus on constructing a general principle that applies to different numerical approaches, and refrain from touching on the issue of regularization. The terms in Eq.( 6) may also be viewed as regularizers, or deformation priors, as they penalize log-Jacobian values that deviate from zero. However, it is well-known in the grid generation field that an integral constraint on the Jacobian alone does not generally guarantee a globally smooth grid [START_REF] Liseikin | Grid Generation Methods[END_REF], so the smoothness of the resulting maps deserves further study. Lastly, to complete our discussions on implementation issues, let us provide the gradient descent direction contributed by the symmetric KL distance term in eq [START_REF] Collins | Automated 3D nonlinear deformation procedure for determination of gross morphometric variability in human brain[END_REF], via its Euler-Lagrange equation. To this end, let us denote Co ij , the matrix cofactor for the (i,j) th component of the Jacobian matrix Dh, we then obtain its Euler-Lagrange equation, using standard Calculus of variations with respect to the i-th coordinate as follows

( ) { } ( ) ( ) . | ) ( | ) ( ) ( , 0 ) ( ) ( | | 1 | ) ( | log 1 1 x Dh x Co x Dh x Co x Dh x Dh x T ij ij j j = = - + ∂ ∂ - ∑ (9)
Figure 1 shows a numerical example demonstrating the proposed approach implemented using a gradient descent projection method. Here, we matched the sequential MRI images (both of size 128 by 128 by 128) from a single subject diagnosed with semantic dementia. The figure showed three 2D slices, plotting the source, target, as well as the deformed source images. As in other TBM approaches, the technique is valuable as the pattern of atrophy is computed automatically, without interactive specification of regions of interest.

Conclusion

This paper is the first effort to systemically examine the relationship of image registration, information theory, and grid generation. While information-theoretic measures such as mutual information [START_REF] Agostino | A viscous fluid model for multimodal non-rigid image registration using mutual information[END_REF][START_REF] Kim | Mutual information for automated unwarping of rat brain autoradiographs[END_REF][START_REF] Lorenzen | Model based symmetric information theoretic large deformation multi-modal image registration[END_REF][START_REF] Twining | A Unified Information-Theoretic Approach to Groupwise Non-Rigid Registration and Model Building[END_REF], f-divergence [START_REF] Pluim | f-information measures in medical image registration[END_REF] and Jensen-Rényi divergence [START_REF] Leow | Asymmetrical Effects of Lithium on Brain Structure Mapped in Healthy Individuals[END_REF] have been popular measures to describe intensity correspondences in nonlinear image registration, it is much less common to appeal to statistical divergence measures in analyzing deformation fields. The proposed formulation calculates the KL-divergence between the deformation and the identity map, treating them as density functions defined on images. Unlike approaches employing conventional continuum mechanics (e.g., Lamé coefficients), our formulation is unbiased and parameter-free (the only parameter involved is the weight, which may be viewed as a Lagrange multiplier). The approach does not therefore make strong assumptions to model the deformation process. In contrast, the commonly employed elasticity theory assumes that displacements or velocity fields to conform to the law of elasticity, which most likely does not accurately describe any real brain deformation process over time.

The proposed framework helps explain the need for log-transformation of Jacobian values in TBM studies, which is ubiquitous and essential in analyzing tissue shrinkage/expansion. Symmetrization is also fundamental to securing inverse-consistency, an important property in image registration. Lastly, this framework is also consistent with a large-deformation approach [START_REF] Avants | Geodesic estimation for large deformation anatomical shape averaging and interpolation[END_REF][START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF][START_REF] Miller | Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms[END_REF], as any one-to-one and onto diffeomorphism remains in the solution space (since its Jacobian is positive and finite everywhere, and thus is theoretically attainable). In the future, we intend to further investigate this approach, which we believe provides a new perspective in nonlinear registration and connects it to other fields of mathematics.

Fig. 1 .

 1 Fig. 1. Image registration using the unbiased deformation algorithm, with the sum of the squared intensity difference (SSD) as the cost function. 3D scans from a semantic dementia patient imaged at two time points were nonlinearly registered to estimate the profile of volumetric change. The later scan (first column) was chosen to be the source image. Progressive brain atrophy is observed in these axial slices.


Cristina Teles

Andrew Smith

Gordon Ramage

Sue Lang
email: sue.lang@gcu.ac.uk.

Identification of clinically relevant viridans group streptococci by phenotypic and genotypic analysis

Keywords: Viridans groups streptococci species identification Streptococcus, MLSA, API, VITEK, RFLP, sodA

Purpose: Two phenotypic and three molecular methods were assessed for their ability to identify viridans group streptococci (VGS) to the species level.

Methods: A panel of 23 clinical isolates, comprising strains isolated from infective endocarditis, blood cultures, pleural and peritoneal fluid, and 19 type/reference strains were analyzed. Identification was performed using two conventional phenotypic methods, API ® rapid ID 32 Strep and VITEK ® 2 system; and genotypic analysis of the nucleotide sequence of the housekeeping gene sodA, restriction patterns generated by restriction fragment length polymorphism (RFLP) of the 16S rRNA gene and multilocus sequence analysis (MLSA) of 7 housekeeping genes.

Results:

The API ® rapid ID 32 Strep accurately speciated 79% of the strains assessed, while the VITEK ® 2 generated a successful identification for 55%, presenting limitations particularly with regard to species belonging to the mitis group. RFLP of the 16S rRNA gene correctly speciated 24% of the strains, having failed to allocate a species for 36% of the isolates examined. In contrast, sequence analysis of the sodA gene provided a correct identification for 95% of the strains assessed, while identification using the MLSA technique was unsuccessful due to practical limitations.

Conclusions:

The results generated herein indicate that no single methodology can be used to provide an accurate identification to the species level of all VGS although, nucleotide sequence analysis of the sodA gene proved to be useful in providing reliable speciation.

INTRODUCTION

Viridans group streptococci (VGS) encompass four phylogenetic clusters: the mitis, mutans, salivarius and anginosus groups, which are part of the natural human flora, being isolated mainly from the oral cavity, gastrointestinal and genitourinary tract [START_REF] Alam | PCR-Based methods for genotyping viridans group streptococci[END_REF]. The pathogenic potential of the species comprising the VGS is variable, with the natural competence of these organisms giving rise to a genetic exchange with Streptococcus pneumoniae, a species with recognized virulence traits [START_REF] Hakenbeck | Molecular Biology of Streptococci[END_REF]. It is believed that this interspecies gene transfer is responsible for the emerging tolerance or resistance to β-lactam antimicrobial agents amongst VGS, particularly S. mitis [START_REF] Amoroso | All detectable high-molecular-mass penicillin-binding proteins are modified in a high-level beta-lactam-resistant clinical isolate of Streptococcus mitis[END_REF]. VGS can behave as opportunistic pathogens causing diseases that range from mild infections to life threatening conditions. S. mutans for instance, has a strong association with the development of dental caries [START_REF] Oho | Simple and rapid detection of Streptococcus mutans and Streptococcus sobrinus in human saliva by polymerase chain reaction[END_REF], while S. sanguinis, S. oralis and S. mitis are responsible for 40-60% of native-valve infective endocarditis [START_REF] Beynon | Infective endocarditis[END_REF][START_REF] Chen | Identification of clinically relevant viridans group streptococci by sequence analysis of the 16S-23S ribosomal DNA spacer region[END_REF]. Members of the anginosus group are causative agents of pyogenic abscesses in a number of body sites [START_REF] Claridge | Streptococcus intermedius, Streptococcus constellatus, and Streptococcus anginosus ("Streptococcus milleri group") are of different clinical importance and are not equally associated with abscess[END_REF]. Other opportunistic infections caused by VGS include septicaemia in neutropenic patients [START_REF] Tunkel | Infections caused by viridans streptococci in patients with neutropenia[END_REF], and occasionally cases of meningitis have been reported [START_REF] Shenep | Viridans-group streptococcal infections in immunocompromised hosts[END_REF].

Currently, no recognised "gold standard" exists that is able to identify many of the species satisfactorily, which has led to an inadequate reporting of streptococcal infections, and potentially has serious implications for diagnosis and treatment [START_REF] Boudewijns | 16S rRNA gene sequencing and the routine clinical microbiology laboratory: a perfect marriage?[END_REF][START_REF] Petti | The role of 16S rRNA gene sequencing in identification of microorganisms misidentified by conventional methods[END_REF][START_REF] Yuan | Clinical impact associated with corrected results in clinical microbiology testing[END_REF].

Accurate identification of VGS to the species level is therefore important to characterise the pathogenic potential of individual species and monitor trends in antimicrobial susceptibility and emerging infections, such as septicaemia in immunocompromised individuals [START_REF] Haanpera | Identification of alpha-hemolytic streptococci by pyrosequencing the 16S rRNA gene and by use of VITEK 2[END_REF].

Although presenting recognized limitations, phenotypic tests such as the API ® rapid ID 32 Strep and the VITEK ® 2 system remain the techniques most frequently used in clinical laboratories, achieving a speciation on the basis of the organisms biochemical properties [START_REF] Jayarao | Comparative evaluation of Vitek gram-positive identification system and API Rapid Strep system for identification of Streptococcus species of bovine origin[END_REF][START_REF] Lupetti | Rapid identification and antimicrobial susceptibility profiling of Gram-positive cocci in blood cultures with the Vitek 2 system[END_REF]. The demand for simple, effective and reliable methods for species identification both in the clinical and research laboratories has resulted in an increasing utilization and development of molecular based techniques [START_REF] Innings | Identification of 43 Streptococcus species by pyrosequencing analysis of the rnpB gene[END_REF]. Nucleotide analysis of the sodA gene, encoding a manganese-dependent superoxide dismutase enzyme (Mn-SOD), has been proposed as an accurate identification tool [START_REF] Hoshino | Use of phylogenetic and phenotypic analyses to identify nonhemolytic streptococci isolated from bacteremic patients[END_REF][START_REF] Poyart | Identification of streptococci to species level by sequencing the gene encoding the manganese-dependent superoxide dismutase[END_REF]. Other genotypic tests previously described include analysis of restriction patterns of the total chromosome (pulsed-field gel electrophoresis) [START_REF] Baseggio | Strain differentiation of isolates of streptococci from bovine mastitis by pulsed-field gel electrophoresis[END_REF][START_REF] Mineyama | DNA fingerprinting of isolates of Streptococcus mutans by pulsed-field gel electrophoresis[END_REF] or within specific genes [START_REF] Sato | Identification of mutans streptococci by restriction fragment length polymorphism analysis of polymerase chain reaction-amplified 16S ribosomal RNA genes[END_REF]. RFLP of ribosomal RNA encoding genes has been used as an identification method for several bacterial and yeast species including Leuconostoc, Malassezia and Streptococcus [START_REF] Jang | A rapid method for identification of typical Leuconostoc species by 16S rDNA PCR-RFLP analysis[END_REF][START_REF] Mirhendi | A simple PCR-RFLP method for identification and differentiation of 11 Malassezia species[END_REF][START_REF] Rudney | Restriction endonuclease-fragment polymorphisms of oral viridans streptococci, compared by conventional and fieldinversion gel electrophoresis[END_REF]. More recently, a MLSA technique has been developed for identification of VGS based on the nucleotide sequence analysis of seven housekeeping genes [START_REF] Bishop | Assigning strains to bacterial species via the internet[END_REF]. Using an electronic database comprising the nucleotide sequences acquired from 420 VGS strains, this methodology reportedly allows unknown strains to be assigned to viridans species via the internet [START_REF] Bishop | Assigning strains to bacterial species via the internet[END_REF].

Here we report the applicability of three molecular based techniques, as well as, two commercially available methods, widely used in routine clinical laboratories, in the identification of VGS to the species level.

MATERIALS AND METHODS

Bacterial strains and culture conditions. A collection of 23 clinical isolates, comprising 14 streptococci recovered from patients diagnosed as Duke definite for infective endocarditis, 7 strains isolated from blood cultures, 1 pleural fluid isolate and 1 peritoneal fluid (Glasgow Caledonian University culture collection, Glasgow Dental School, Barts and The London NHS Trust) were examined (Table 1). For assessment of the accuracy of the identification provided by each method 10 type and 9 reference strains were included (Table 1). The strains designated ATCC and NCTC were obtained from the American Type Culture Collection (Manassas, USA) and National Collection of Type Cultures (London, UK), respectively.

All strains were grown at 37ºC in an atmosphere enriched with 5% CO 2 on Columbia base agar (Oxoid, Hampshire, UK) supplemented with 5% (v/v) defibrinated horse blood (Oxoid) and cultivated in brain heart infusion broth (BHI) (Oxoid).

Phenotypic identification. The API ® rapid ID 32 Strep and VITEK ® 2 system kits were purchased from bioMérieux, UK. Phenotypic analysis by both methods was performed according to the manufacturer's instructions. The API results were obtained as a numerical profile after visual reading and the identification was achieved by manually entering the code into the identification software (apiweb -API ® rapid ID 32 Strep V3.0, bioMérieux, France). Identification through the VITEK ® 2 system was provided by the associated computer software (V1.02), through comparison of the biochemical profile of the strains tested against the biochemical profiles of the strains present in the software database.

Genotypic identification

DNA manipulations. Streptococcal DNA was extracted using the DNeasy ® Tissue kit (Qiagen, West Sussex, UK) and quantified using a ND-1000 spectrophotometer (NanoDrop Technologies Inc, Wilmington, USA).

Polymerase chain reaction (PCR) analysis. DNA amplifications were performed using an Applied Biosystems 2720 thermal cycler (Applied Biosystems, Foster City, USA). All oligonucleotide primer sequences are depicted in table 2. The PCR reactions consisted of 1× NH 4 amplification buffer, 1.5mM MgCl 2 , 1mM dNTP mix, 1µM of each primer, 10ng of genomic DNA and 0.025U/µl Taq polymerase (Bioline, London, UK) in a final volume of 50µl. Amplification of the sodA gene was performed on all strains using an annealing temperature of 51.4ºC.

S. oralis NCTC 11427 T and S. oralis 12601 were selected to assess the MLSA technique for VGS species identification. Genes map, pfl and tuf were amplified with 55ºC annealing temperature, while primers for genes ppaC, pyK and rpoB were annealed at 50ºC, as previously published [START_REF] Bishop | Assigning strains to bacterial species via the internet[END_REF]. PCR amplifications were performed with 30 cycles of denaturation (94ºC, 30 seconds), annealing (XºC, 30 seconds) and extension (72ºC, 30 seconds), followed by a final extension step (72ºC, 5 minutes) [START_REF] Bishop | Assigning strains to bacterial species via the internet[END_REF].

All generated amplicons were resolved by electrophoresis using a 1.0% (w/v) SeaPlaque ® agarose gel (Cambrex, Bio Science Rockland, Inc., USA), excised from the gel and purified using the QIAquick ® gel extraction kit (Qiagen).

The amplification of the 16S rRNA gene was performed with 35 cycles of denaturation (94ºC for 1 minute), annealing (60ºC for 1 minute) and extension (72ºC for 1.5 minutes), followed by a final extension step at 72ºC for 10 minutes [START_REF] Sato | Identification of mutans streptococci by restriction fragment length polymorphism analysis of polymerase chain reaction-amplified 16S ribosomal RNA genes[END_REF]. The resulting products were resolved by electrophoresis using a 1% (w/v) agarose gel (Invitrogen, Paisley, UK), after which the amplicons were recovered and purified from the PCR reaction mixture using QIAquick® PCR Purification Kit (Qiagen).

Sequencing analysis. Nucleotide sequence analysis of the sodA gene was performed based on the nucleotide sequence analysis of a single DNA strand of all VGS examined. Evaluation of MLSA was performed by sequencing both DNA strands of the amplified fragments of the sodA, map, pfl, tuf, ppaC, pyK and rpoB genes of S. oralis NCTC 11427 T and S. oralis 12601 (infective endocarditis isolate). Sequencing was performed using 15ng of each amplicon, 3.5µM primer and 50ºC annealing temperature or 64µM RFLP of the 16S rRNA gene. RFLP of the amplified 16S rRNA genes was carried out as described in previous studies [START_REF] Sato | Identification of mutans streptococci by restriction fragment length polymorphism analysis of polymerase chain reaction-amplified 16S ribosomal RNA genes[END_REF], with modifications regarding the restriction enzyme digestion conditions. Double restriction enzyme digestions were performed with HaeIII and HpaII (New England Biolabs, Hitchin, UK) according to the manufacturer's instructions and the resulting restriction fragments were resolved by electrophoresis using a 2% (w/v) agarose gel (Invitrogen). A library of the restriction patterns was produced for the reference strains and identification of the clinical isolates was achieved using Dice similarity coefficient. A dendrogram was then constructed showing the clusters generated by the Unweighted pair group method with arithmetic mean (UPGMA, BioNumerics software, Applied Maths NV, Belgium).

RESULTS

Species identification through API ® rapid ID 32 Strep system. The API ® rapid ID 32

Strep system provided an accurate identification to the species level for 33 strains (79%), without the need to perform additional tests. Eight strains (19%) were given a correct identification to the group level, these included 7 strains from the mitis group, which presented a number of unusual reactions, and type strain S. sobrinus ATCC 33478 T was given an inconclusive identification, stating sobrinus/S. downei. Since the biochemical range of this method is unable to differentiate between these two species further tests were recommended.

Species identification through the Vitek ® 2 system. The results generated by the Vitek ® 2 accurately speciated 23 strains (55%) while, 9 isolates of the mitis group were allocated to the correct group (21%). The result of all S. oralis and S. mitis was presented as S. mitis/S. oralis with a high level of confidence, since this methodology does not differentiate between these species. Assignment of S. sanguinis reference and type strains to S. mitis/S. oralis occurred in two instances (S. sanguinis NCTC 7863 and S. sanguinis ATCC10556 T ), possibly due to a negative reaction in both arginine tests, which is characteristic of S. mitis and S. oralis strains [START_REF] Facklam | What happened to the streptococci: overview of taxonomic and nomenclature changes[END_REF]. S. parasanguinis ATCC 15912 T was speciated as S. gordonii, presenting a single obvious atypical reaction (methyl-β-D-glucopyranoside positive). No identification was generated for 8 strains (19%) represented in the software database, including isolates from the mitis, salivarius and anginosus group, with 2 (5%) isolates having been misidentified. Species identification through sequence analysis of the sodA gene. The PCR amplification of an internal fragment of the sodA gene was successful for all strains, resulting in a single amplicon with the expected size of approximately 440bp. Sequence analysis of the sodA gene resulted in the correct speciation of 40 strains (95%), with S. sobrinus ATCC 33478 T having been incorrectly assigned to S. thermophilus with 89% nucleotide similarity and failure to generate sequence data occurred for one isolate, despite repeated attempts.

Species identification through RFLP of the 16S rRNA gene. Amplification of the 16S rRNA gene of all VGS strains generated a single amplicon with the expected size of approximately 1500bp. With the similarity threshold set at 90%, analysis of the RFLP revealed 11 clusters incorporating 27 strains (64%), 24% of which were correctly speciated (10/42) and 21% of isolates (9/42) were allocated to the correct VGS group while 19% (8/42) were misidentified (fig. 1). The remaining 15 strains (36%) examined, including species of the mutans, mitis and anginosus group, were excluded from the generated clusters despite the presence of other members of the same species and/or VGS group (fig. 1).

Species identification through MLSA. The PCR amplification of the 7 housekeeping genes was successful for both strains examined, generating single amplicons of the expected sizes. Using the forward primers the nucleotide sequence generated for 3 genes (3/7) of the type strain S. oralis NCTC 11427 T and 2 genes (2/7) of the clinical isolate S. oralis 12601 did not reach the length required by the MLSA software and both strains failed to yield sequence results using the ppaC-F primer. When analyzing the complementary DNA strand using the reverse primers, the ppaC gene of S. oralis NCTC 11427 T and S. oralis 12601 remained below the number of base pairs necessary to proceed with the analysis. In addition, no nucleotide sequence was obtained for the pfl and pyk genes of the type strain and map, pyk, rpoB and sodA genes of the clinical isolate. When modifications were introduced to the published protocol, an improvement was seen for 2 of the genes of interest however, this was not beneficial overall as it caused no improvement on either strand of the ppaC gene and decreased the number of base pairs acquired for the sodA gene. Since the required read length of the 7 housekeeping genes was not obtained for either of the strains examined, no analysis could be performed.

Identification concordance between the methods assessed.

Comparison of the identification allocated to each strain by the techniques assessed, revealed a species level agreement between the two phenotypic methods for 18 strains (43%). In contrast, the results generated by the genotypic methods were concordant for 10 isolates (24%).

The species allocated by the VITEK ® 2 and sequence analysis of the sodA gene were similar for 21 strains (50%) however, when the identifications produced by the VITEK ® 2 were compared with those of the RFLP of the 16S rRNA gene, only 6 strains (14%) were allocated to the same species. The API ® rapid ID 32 Strep and the RFLP of the 16S rRNA gene were concordant for 8 strains (19%) while, the highest levels of agreement were obtained when comparing the API ® rapid ID 32 Strep and sequence analysis of the sodA gene, with 32 isolates (76%) being similarly speciated.

DISCUSSION

The results presented in this study reflect the complex taxonomy of the VGS which has lead to difficulties in defining a "gold standard" technique that would provide accurate identification to the species level. The API ® rapid ID 32 Strep system remains one of the most frequently used methodologies in the clinical laboratory for streptococci identification to the species level [START_REF] Jayarao | Comparative evaluation of Vitek gram-positive identification system and API Rapid Strep system for identification of Streptococcus species of bovine origin[END_REF]. Previous reports demonstrated the low sensitivity of this technique for speciation of members of the anginosus and mitis group [START_REF] Chen | Identification of clinically relevant viridans group streptococci by sequence analysis of the 16S-23S ribosomal DNA spacer region[END_REF][START_REF] Freney | Description and evaluation of the semiautomated 4-hour rapid ID 32 Strep method for identification of streptococci and members of related genera[END_REF][START_REF] Jacobs | Rapid species identification of "Streptococcus milleri" strains by line blot hybridization: identification of a distinct 16S rRNA population closely related to Streptococcus constellatus[END_REF][START_REF] Limia | Comparison of three methods for identification of Streptococcus milleri group isolates to species level[END_REF],

which is in agreement with the results found in the present study as all the discrepant speciations by this method occurred with members of the mitis group. Further limitations of this system, and others of this nature, are associated with the subjectivity inherent with the interpretation of the colorimetric reactions which can interfere with the identification obtained. In addition, the availability of an up-to-date database incorporating the taxonomic re-evaluations of this diverse group of bacteria is crucial in acquiring accurate results [START_REF] Bosshard | Comparison of conventional and molecular methods for identification of aerobic catalase-negative gram-positive cocci in the clinical laboratory[END_REF]. Despite being a simple and non time consuming method, the API ® rapid ID 32 Strep only provides reliable identification to the group level.

The VITEK ® 2 system is an automated machine that is capable of identifying numerous bacteria, including VGS, based on kinetic analysis of the fluorescence, turbidity and colorimetric readings of 43 tests measuring enzymatic activities, carbon source and antibiotic resistance [START_REF] Ligozzi | Evaluation of the VITEK 2 system for identification and antimicrobial susceptibility testing of medically relevant gram-positive cocci[END_REF]. Difficulties in the speciation of VGS members using the VITEK ® 2 have been reported [START_REF] Innings | Identification of 43 Streptococcus species by pyrosequencing analysis of the rnpB gene[END_REF][START_REF] Brigante | Use of the Phoenix automated system for identification of Streptococcus and Enterococcus spp[END_REF] and the results generated by this technique in the present assessment reinforce these findings in that only 55% of the strains examined were accurately speciated. In addition, the supplementary tests required for two of the most clinically relevant VGS species, S. mitis and S. oralis, proves to be a major limitation also extending the time required to obtain identification.

In agreement with previous reports, the results provided by the VITEK ® 2 suggest this methodology may be better suited for identification to the species level of organisms other than VGS, such as β-haemolytic streptococci [START_REF] Haanpera | Identification of alpha-hemolytic streptococci by pyrosequencing the 16S rRNA gene and by use of VITEK 2[END_REF].

The sequence analysis of an internal fragment of the sodA gene resulted in an accurate speciation of 95% of all strains assessed, revealing the high reliability of this method for VGS speciation. Traditionally, the sequencing of this housekeeping gene has been carried out by cloning the amplified fragment into a vector prior to sequencing [START_REF] Hoshino | Use of phylogenetic and phenotypic analyses to identify nonhemolytic streptococci isolated from bacteremic patients[END_REF][START_REF] Poyart | Identification of streptococci to species level by sequencing the gene encoding the manganese-dependent superoxide dismutase[END_REF][START_REF] Poyart | Characterization of superoxide dismutase genes from gram-positive bacteria by polymerase chain reaction using degenerate primers[END_REF]. In this study, elimination of this step made this technique less time-consuming and more cost-effective, without significant consequence to its efficiency. GenBank includes a comprehensive and freely available microbial genomic database which, despite not guaranteeing the accuracy and quality of the data deposited, is an up-to-date database which incorporates the sodA nucleotide sequences of VGS reference/type strains and clinical isolates [START_REF] Hoshino | Use of phylogenetic and phenotypic analyses to identify nonhemolytic streptococci isolated from bacteremic patients[END_REF][START_REF] Poyart | Identification of streptococci to species level by sequencing the gene encoding the manganese-dependent superoxide dismutase[END_REF][START_REF] Poyart | Characterization of superoxide dismutase genes from gram-positive bacteria by polymerase chain reaction using degenerate primers[END_REF]. These factors have contributed to the high levels of reliable identification generated by this methodology.

The sequence analysis of the 16S rRNA gene as been widely applied for speciation of streptococci having greatly contributed to a better understanding of the classification of such a complex bacterial group [START_REF] Kawamura | Determination of 16S rRNA sequences of Streptococcus mitis and Streptococcus gordonii and phylogenetic relationships among members of the genus Streptococcus[END_REF]. The high level of nucleotide similarities between closely related species, particularly S. oralis, S. mitis and S. pneumoniae soon became evident, questioning the discriminatory capability of this methodology [START_REF] Hoshino | Use of phylogenetic and phenotypic analyses to identify nonhemolytic streptococci isolated from bacteremic patients[END_REF]. Restriction enzyme digestion of the 16S rRNA gene has been suggested as an alternative method to gene sequence analysis for streptococcal speciation, potentially offering more precise results [START_REF] Sato | Identification of mutans streptococci by restriction fragment length polymorphism analysis of polymerase chain reaction-amplified 16S ribosomal RNA genes[END_REF]. The use of a single restriction enzyme (HpaII) for discrimination of mutans group species was previously described [START_REF] Sato | Identification of mutans streptococci by restriction fragment length polymorphism analysis of polymerase chain reaction-amplified 16S ribosomal RNA genes[END_REF] however, in that study, not all strains presented distinctive patterns requiring a second restriction endonuclease (HaeIII) to allow differentiation. Given the number and diversity of strains being assessed in the present evaluation, a simultaneous double restriction reaction was applied. Traditionally, species identification was performed by visual comparison of the restriction patterns resolved by agarose gel electrophoresis however, computer analysis is less subjective and a more practical approach for high throughput of strains [START_REF] Rudney | Restriction endonuclease-fragment polymorphisms of oral viridans streptococci, compared by conventional and fieldinversion gel electrophoresis[END_REF]. Analysis of the constructed dendrogram demonstrated the limited capacity of RFLP of the 16S rRNA gene to speciate VGS (24%) and therefore, low concordance with the other methodologies evaluated. This unsatisfactory results may have been influenced by the low resolution, or inability to visually detect, small restriction fragments. Other limitations of this technique have also been reported. These included the inability to differentiate between S. parasanguinis and S. mitis, despite the use of three separate restriction endonuclease combinations, due to the sequence homology of the 16S rRNA gene shared amongst strains within a VGS group [START_REF] Rudney | Restriction endonuclease-fragment polymorphisms of oral viridans streptococci, compared by conventional and fieldinversion gel electrophoresis[END_REF].

Based on the findings generated using the panel of strains assessed, RFLP of the 16S rRNA genes was considered unsuitable for streptococcal identification to the species level.

Unlike the MLST technique successfully applied in the genetic differentiation between serotypes of S. aureus and S. pneumoniae, MLSA has recently been described for VGS species identification based on the sequence data acquired for 7 housekeeping genes [START_REF] Bishop | Assigning strains to bacterial species via the internet[END_REF]. The results obtained in the present study suggest that MLSA of VGS species, at present is not a technique sufficiently robust to allow inter-lab variability and cannot be used routinely.

In conclusion, no individual method assessed was able to unequivocally assign the correct species to all VGS strains. The difficulties in achieving accurate speciation were particularly noticeable in members of the mitis group, which correlates with the genetic homology shared among the species included in this VGS group. Comparison of the data provided indicated that nucleotide sequence analysis of the sodA gene is a useful and reliable technique for VGS species identification.

  primer and 45ºC annealing temperature by The Sequencing Service (College of Life Sciences, University of Dundee, UK) using Applied Biosystems Big-Dye Ver 3.1 chemistry on an Applied Biosystems model 3730 automated capillary DNA sequencer. Species identification was achieved through nucleotide comparison against published sequences available in the microbial genomes database using the Basic Local Alignment Search Tool (BLAST) available online (http://www.ncbi.nlm.nih.gov/BLAST/).

Table 2 -

 2 PCR primers to amplify sodA, 16S rRNA, map, pfl, tuf, ppaC, pyK and rpoB genes. The underlined degenerate base sites correspond to positions within the gene that may differ between species.

				Expected
			Gene	product
	Primers	Sequence (5' to 3')		Origin
			amplified	size
				(bp)
	sodA-F TRCAYCATGAYAARCACCAT *
			sodA	436	[18]
	sodA-R ARRTARTAMGCRTGYTCCCARACRTC *
	8UA	AGAGTTTGATCCTGGCTCAG	16S
				1521	[22]
	1492R	TACGGGTACCTTGTTACGACTT	rRNA
	map-F	GCWGACTCWTGTTGGGCWTATGC
			map	459
	map-R	TTARTAAGTTCYTTCTTCDCCTTG
	pfl-F	AACGTTGCTTACTCTAAACAAACTGG
			pfl	458
	pfl-R	ACTTCRTGGAAGACACGTTGWGTC*
	tuf-F	GTTGAAATGGAAATCCGTGACC
			tuf	531
	tuf-R	GTTGAAGAATGGAGTGTGACG
				[26]
	ppaC-F GACCAYAATGAATTYCARCAATC*
			ppaC	683
	ppaC-R TGAGGNACMACTTGTTTSTTACG*
	pyk-F	GCGGTWGAAWTCCGTGGTG*
			pyK	557
	pyk-R	GCAAGWGCTGGGAAAGGAAT*
	rpoB-F AARYTIGGMCCTGAAGAAAT*
			rpoB	740
	rpoB-R TGIARTTTRTCATCAACCATGTG*

*

evaluation of the methods accuracy.