Jérôme Gouin

Catherine Lerouge

Yves Deschamps

David Widory

Fatima Laggoun-Défarge

Mariana Yossifova

Yotzo Yanev

Petrological and geochemical characterization of Ge-bearing coals from the eastern Rhodopes, Bulgaria

Keywords: coal geochemistry, germanium, sulphur isotopes, eastern Rhodopes, Bulgaria

The coals from Pchelarovo and Vulce Pole are perhydrogenated coals showing high sulphur content (3.2 -6.2 wt%) and significant Ga and Ge contents (0 -480 ppm) in organic matter (OM). The coals from Medenbuk could be perhydrogenated coals modified by weathering processes. They show low sulphur content (1.5 wt%), high homogeneous Ge contents in OM (~2500 ppm) and an unusual enrichment in Ge (up to 0.8 wt%) and V (2.8 wt%) in weathered zones of framboïdal pyrite. Sulphur isotopic study of pyrite and organic sulphur provides evidence of bacterial reduction of dominantly seawater sulphates. However taking account the geological context of the three basins, the introduction of volcanic fluids in the basins of Pchelarovo and Vulce Pole during their formation cannot be excluded.

INTRODUCTION

Coal is the main energy resource in Bulgaria and is the subject of numerous investigations on metal contents on the primary resource and on ashes after combustion in power stations. Coals from Pchelarovo, Vulce Pole and Medenbuk were previouly studied for their Ge content [START_REF] Eskenazy | Geochemistry of germanium in Bulgarian coal[END_REF]. Investigations on sulfur content and distribution, and on relationships between organic matter and metals are good indicators of the geological deposition. This abstract presents preliminary petrological and geochemical data of coals from the three coal deposits, in the aim to understand their conditions of deposition.

GEOLOGICAL SETTING

The Pchelarovo, Vulce Pole and Medenbuk coal deposits are located in the Eastern Rhodopes, resulting from the continental Macedonian -Rhodope -North-Aegean magmatic arc [START_REF] Harkovska | Paleogene magmatism in the Central Rhodopes area, Bulgaria -A review and new data[END_REF], formed after the Late Cretaceous-Palaeocene collision of the Serbo-macedonian and Rhodopes Massifs with the Pelagonian microplate. Coal-bearing sediments were deposited within intra-orogenic ex-tensional basins during the final stage of continental collision in the Paleogene.

Pchelarovo

The Pchelarovo coal deposit is located about 20 km north of Kardzhali. Coal-sediments belong to the Pchelarovo depression quite close to the southeast of the Borovitsa volcanic area. The basement of the basin is composed of immature conglomerates and sandstones, reworking the metamorphic basement (Central Rhodopean dome and metamorphosed Palaeozoic basement). The basin is infilled by Paleogene sediments composed of pyroclastic and volcanic rocks of upper Eocene and Oligocene ages. The petrology of the sandstones intercalated with coal layers at Pchelarovo shows a high proportion of devitrified hyaloclastic components in the sediments, indicating that coal deposition is contemporaneous with explosive volcanism whose volcanic centers were close to the northeastern border of the Borovitsa caldera. An intercalated marl-calcareous layer with microfossils indicates a marine transgression. A previous study of advanced argillic alteration (AAA) in Bulgaria has shown that volcanic rocks of the area are strongly hydrothermally altered (Kunov, 1999). Complementary isotopic (S, O, H) data on alunite indicate that AAA is of magmatic-hydrothermal origin with an important participation of external seawater-derived fluids [START_REF] Lerouge | Constraints of stable isotopes on the origin of alunite from advanced argillic alteration systems in Bulgaria[END_REF]. Coals occur as 5cm-thick vitrinite layers rich in pyrite (samples 6 to 9), or as thin seams (11)(12). The sample 7 is characterized by small perpendicular fractures infilled by pyrite.

Vulce Pole

The coal samples were collected close to the Vulce Pole village in a tectonic depression 25 km long in east-west direction and 3-4 to 10 km wide, filled by Oligocene sediments. Coals occur as 10 to 50 cm-thick seams (samples 29 and 31). The basement of the basin consists essentially of gneisses and amphibolites. Coalbearing sediments of the basin are dominantly conglomerates and sandstones with intercalated volcanic tuffs and rare limestone at the bottom of the series, and clay layers at the top. The petrology of conglomerates and sandstones indicates immature sediments with a dominant contribution from the reworked basement (gneiss, granite) and an important component of devitrified volcanic clasts, suggesting that coal deposition is contemporaneous with explosive volcanic episode. Matrix of the elements is not abundant; carbonates are absent or rare, nodules of pyrite infill partially the porosity of some layers of sandstone and conglomerates suggesting sedimentation rather in a continental domain with late circulation of H 2 S-bearing fluids. The volcanic centre close to the basin is the Madjarovo complex (120 km 2 ) in the south, dated by the 40 Ar/ 39 Ar method on biotite and sanidine at around 32-33 Ma (Singer and Marchev, 2002). It induced intense magmatichydrothermal activity, marked by minor barren advanced argillic alteration and a welldeveloped epithermal adularia-sericite alteration associated with a base-metal-Au-Ag mineralization (Singer and Marchev, 2002).

Medenbuk

The coal sample (58) was collected from an outcrop of Paleogene sediments, located 20 km from Medenbuk. The basement of the basin is dominantly represented by gneisses and amphibolites of Proterozoic age. The coal-bearing sediments are close to the border and the basement of the basin; they are immature sandstones reworking exclusively components from the basement, consisting of gneisses, bioclastrich carbonates but also mafic rocks marked by presence of epidote, amphibole and clinopyroxene. The matrix is abundant and represented by micritic carbonates, strongly suggesting deposition under relatively deep marine conditions. The coal was present as small lenses, then partly reworked in a fracture. Nodules of pyrite are observed in sandstones directly above coal lenses.

RESULTS

Coal characterization

HI, OI, Tmax, S1 and S2 were analysed by Rock-Eval® RE6 pyrolysis according to standard methods. Vitrinite reflectance measurements were performed using a Leica DMR XP microscope under reflected light, following ICCP procedures. The coals from Pchelarovo are composed of dominantly collotelinite with minor humotelinite and semifusinite, whereas the coals from Vulce Pole and Medenbuk are composed of dominantly humotelinite with some collotelinite. The vitrinite from the three basins shows high carbon contents (79-86 wt%), consistent with a bituminous high to medium volatile coal rank. However the vitrinite reflectance and Tmax values are low, rather indicating a subbituminous rank (Table 1). The coals from Vulce Pole and Pchelarovo are characterized by a bright green-yellow fluorescence and by abnormally high HI values and relatively low OI. These characteristics are typical of perhydrogenated coals (Jimenez et al., 1998). The sample 9 from Pchelarovo has OI and HI values characteristic of humic coals (Fig. 1). The coal from Medenbuk is very different with a low fluorescence intensity high OI and low HI values. These characteristics may be interpreted by two ways. First, the sample is an immature humic coal. Second, it could be a perhydrogenous coal comparable to that from Vulce Pole and Pchelarovo, which could have been altered through hydrothermal processes or weathering. The presence of secondary gypsum and iron oxi-hydroxides in the coal lenses and the low fluorescence are in favor of the second assumption.

S, Ga, Ge and V distribution

Organic matter from the three coal deposits was analyzed for C, S, Fe and metal contents, using a Camebax SX50 electron microprobe. S contents are high whereas Fe contents are low, showing that sulphur can be attributed to organic matter. Organic sulphur contents show that Medenbuk (1.5 wt%) is a medium-sulphur coal and Pchelarovo (3.2 -6.2 wt%) and Vulce Pole (4.3 wt%) are high-sulphur coals. They may be compared with coals overlain by a marine roof in which most of the sulphur is derived from seawater sulfate [START_REF] Chou | Geochemistry of sulfur in coal[END_REF]. Ge and Ga contents were measured in vitrinite of the three ore deposits, with a 20 kV acceleration voltage and a current beam of 20nA. The ranges of Ga content in the three ore deposits are quite similar: 0-260 ppm at Pcelarovo and Vulce Pole, and 0-330 ppm at Medenbuk. The ranges of Ge content at Vulce Pole and Pchelarovo are similar (0-480 ppm). On the other hand, Ge contents in OM from Medenbuk are high (2418 ± 1033 ppm, n = 12) and associated with significant V contents (0-1100 ppm). Unusually high Ge and V contents (respectively up to 0.8 and 2.8 wt%) were measured in iron oxides (Fig. 2) resulting of the alteration of primary framboidal pyrites.

Sulphur isotopes

Sulphur isotope analysis of pyrite, organic matter and secondary gypsum was performed using a CFIRMS system coupled with Flash EA and gas bench. Results are given in Table 2.

Table 2: Sulphur isotopic data on organic matter and sulphides from the three coal deposits.

Pchelarovo

Organic sulphur and syn-sedimentary pyrite have a relatively narrow range of δ 34 S. Independent of the geological context, these values may be interpreted by two ways: 1) bacterial reduction of marine sulphates in a relatively closed system, and 2) input of magmatichydrothermal fluids in the basin during its formation. Taking into account the contemporaneous volcanic activity of the basin (δ 34 S fluid ~ 0-3 ‰, [START_REF] Lerouge | Constraints of stable isotopes on the origin of alunite from advanced argillic alteration systems in Bulgaria[END_REF] and the presence of marine sediments at Pchelarovo, the sulphur origin could be related to both circulation of volcanic-magmatic and seawater-derived fluids during the deposition of the sediments, with the basin evolving as a relatively closed system. The δ 34 S of secant sulphides crosscutting the stratification is around +3 ‰ and is not significantly different from the δ 34 S of the synsedimentary sulphides, suggesting either a partial remobilization of the syn-sedimentary sulphides, or a persistence of circulation of the same type of fluids during sedimentation and diagenesis.

Vulce pole

The δ 34 S of the pyrite infilling the porosity of sandstones varies between +5.6 and +10.2 ‰ (CDT), whereas δ 34 S organic matter is close to +2.2 ‰. Values are heterogeneous and higher than the δ 34 S of the magmatic fluids associated with the Madjarovo volcanic activity. The data cannot exclude the participation of magmatic fluids in the basin, but they preferentially indicate bacterial reduction of sulphates present in the fluids of the basin. The positive and heterogeneous values could result from evolution of sulphate-reducing bacteria in isolated microenvironments.

Medenbuk

The δ 34 S of syn-sedimentary pyrite in coal can be deduced from the secondary gypsum produced by weathering of pyrite. The δ 34 S value of gypsum is -28.5 ‰; that is quite close to the δ 34 S of the pyrite infilling the porosity of sandstone. The δ 34 S of vitrinite is slightly higher than δ 34 S of pyrite, as usual at Pchelarovo and most of coal deposits. This homogeneity of the δ 34 S and the highly negative values are consistent with sedimentary sulphides resulting from the bacterial reduction of marinederived sulphates in an open system for sulphates, and also in good agreement with petrological features of the studied samples.

DISCUSSION

The three coal basins were deposited during the late intra-orogenic extensional phase of the continental collision in the Paleogene. Coals are rich in sulphur incorporated in organic matter and sulphides, and in metals. Coal from Medenbuk shows particularly high Ge content. Sulphur isotopes in organic matter and pyrite are used to interpret fluid origins. The δ 34 S values of organic sulphur and pyrite range from -30 to +10 ‰ (CDT), consistent with sulphate bacterial reduction. However the complexity of the isotopic data requires to take into account the geological context. The three coal basins are different:

1) Depositional environment: The Medenbuk and Pchelarovo basins contain marine sediments just above coal seams. A part of sulphur (in organic matter and sulphides) results from the bacterial reduction of seawaterderived sulphates, in an open system at Medenbuk and in a closed system at Pchelarovo.

2) Contemporaneous volcanism: Pchelarovo and Vulce pole are close to volcanic centres, Medenbuk is far from any volcanism. Positive δ 34 S values may be partly attributed to volcanic fluids involved in the basin at Pchelarovo and Vulce Pole.

Figure 1 :

 1 Figure 1: Location of Bulgarian coal samples in Van Krevelen's HI-OI diagram.

Figure 2 :

 2 Figure 2: Backscattered electron image of partially oxidized framboidal pyrites showing high Ge and V contents.

Table 1 : Rock-Eval pyrolysis, vitrinite reflectance, and C-S data measured by EPMA with a 10 kV accel- eration voltage and a current beam of 20 nA.

 1

	N° C	S	S2 (mg HC	Tmax	Ro	Mean
	%	%	/ g rock)	(°C)	(%)	Ro (%)
			Pchelarovo
	6 84 4.5	75.75	419	0.41
	9 80 6.2	30.61	415	0.47	0.46
	11 80 3.8	90.46	413	0.45
	12 86 3.2	93.46	414	0.49
			Vulce Pole
	29 83 4.3	141.33	370	0.27	0.30
	31 84 4.3	163.88	374	0.33
			Medenbuk
	58 79 1.5	32.01	404	0.32

ACKNOWLEDGEMENTS

This work was financially supported by the BRGM research projects and by the ECO-NET program (10161XM) between the BRGM (France), the Geological Institute of Sofia, Bulgarian Academy of Sciences (Bulgaria) and the University of Skopje (Macedonia).


Maxime Crochemore
email: maxime.crochemore@kcl.ac.uk

Roberto Grossi
email: grossi@di.unipi.it

Juha Kärkkäinen
email: juha.karkkainen@cs.helsinki.fi

Gad M Landau
email: landau@cs.haifa.ac.il

Computing the Burrows-Wheeler transform in place and in small space ✩

We introduce the problem of computing the Burrows-Wheeler Transform (BWT) using small additional space. Our in-place algorithm does not need the explicit storage for the suffix sort array and the output array , as typically required in previous work . It relies on the combinatorial properties of the BWT, and runs in O(n 2 ) time in the comparison model using O(1) extra memory cells, apart from the array of n cells storing the n characters of the input text. We then discuss the time-space trade-off when O(k •σ k ) extra memory cells are allowed with σ k distinct characters, providing an O((n 2 /k + n) logk)-time algorithm to obtain (and invert) the BWT. For example in real systems where the alphabet size is a constant, for any arbitrarily small > 0, the BWT of a text of n bytes can be computed in O(n -1 logn) time using just n extra bytes.

Introduction

The Burrows-Wheeler Transform [START_REF] Burrows | A block-sorting lossless data compression algorithm[END_REF] (known as BWT) of a text string is at the heart of the bzip2 family of text compressors, and finds also applications in text indexing and sequence processing. Consider an input text string T ≡ T [0 . . n -1] and the set of its suffixes T i ≡ T [i . . n -1] (0 ≤ i < n) under the lexicographic order, where T [n -1] is an endmarker character $ smaller than any other character in T . The alphabet from which the characters in T are drawn can be unbounded.

A classical way to define the BWT uses the n circular shifts of the text T = mississippi$ as shown in the first column of Table 1. We perform a lexicographic sort of these shifts, as shown in the second column: if we mark the last character from each of the circular shifts in this order, we obtain a sequence L of n characters that is called the BWT of T . Its relation ✩ A preliminary version of the results in this paper appeared in [START_REF] Crochemore | A constant-space comparison-based algorithm for computing the Burrows-Wheeler transform[END_REF]. The work of the second author has been supported in part by the Italian Ministry of Education, University, and Research (MIUR) under PRIN 2012C4E3KT project. The work of the third author has been supported by the Academy of Finland Grant 118653 (ALGODAN). The work of the fourth author has been partially supported by the National Science Foundation Award 0904246, Israel Science Foundation Grants 347/09 and 571/14, Yahoo, Grant No. 2008217 from the United States-Israel Binational Science Foundation (BSF) and DFG. to suffix sort is well known, as illustrated in the third column: the rth character in L is T [ j -1] if and only if T j is the rth suffix in the sort (except the borderline case j = 0, for which we take T [n -1] as character).

As it can be seen in the example of Table 1, the BWT produces a text of the same length as the input text T . The transform is reversible since it is a one-to-one function when the input text is terminated by an endmarker $. Thus, not only we can recover T from L alone, but typically L is more compressible than T itself using 0th-order compressors [START_REF] Manzini | An analysis of the Burrows-Wheeler transform[END_REF].

There are now efficient methods that convert T to L and vice versa, taking O (n log n) time for unbounded alphabets in the worst case [START_REF] Adjeroh | The Burrows-Wheeler Transform: Data Compression, Suffix Arrays, and Pattern Matching[END_REF]. 1 The BWT is also a key element of some compressed text indexing implementations due to the small amount of space it requires: some examples are the solution by Ferragina and Manzini [START_REF] Ferragina | Indexing compressed text[END_REF] or that by Grossi et al. [START_REF] Grossi | High-order entropy-compressed text indexes[END_REF], where the transform is associated with the techniques of wavelet trees and of succinct data structures using rank-select queries on binary sequences [START_REF] Munro | Tables[END_REF].

One of the prominent applications of the BWT is for software dealing with Next Generation Sequencing, where millions of short strings, called reads, are mapped onto a reference genome. Typical and popular software of this type are Bowtie [START_REF] Langmead | Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[END_REF], BWA [START_REF] Li | Fast and accurate long-read alignment with Burrows-Wheeler transform[END_REF] and SOAP2 [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF]. Here it is crucial that the genome is indexed in a compact manner to get reasonable running time. Space issues for computing the BWT are thus relevant: frequently the input data is so large that the input text T stays in main memory while any additional data structure of similar size cannot fit in the rest of the main memory [START_REF] Wing-Kai Hon | Breaking a time-and-space barrier in constructing full-text indices[END_REF].

All the previous work for computing the BWT of T relies on the fact that (a) we need first to store the suffix sorting of T (also known as suffix array [START_REF] Manber | Suffix arrays: a new method for on-line string searches[END_REF]), thus occupying n memory cells for storing integers, and (b) we need to output the BWT in another array storing n characters. Motivated by these observations, we want to study the case in which (a) and (b) are avoided, thus saving on the space occupied by them.

In this paper, our goal is to obtain the BWT by directly permuting T and using just O (1) memory cells, i.e., we aim at an in-place algorithm for computing the BWT. We consider the model in which the text T is stored as an array of n entries, where each entry stores exactly one character of T . Note that storing an integer usually takes more space than a character, so we assume that only the characters of T can be kept in the array T . Moreover, T is not read-only but it can be modified at any time, and just O (1) additional memory cells (besides T ) can be kept for storing auxiliary information. 2Note that our model represents some realistic situations in which one has to handle large text collections, or large genomic sequences, without relying on extra memory for (a) and (b). Hence it is crucial to maximize the amount of data that can fit into main memory: not storing explicitly (a) and (b) permits to save space, which is typically regarded as taking more than half of the total space required. For instance, DNA sequences are stored by using 2 bits per character and machine integers take 64 bits. Here we just need 2n bits to store the (genomic) text and save the 64n bits used for storing the intermediate suffix sorting in (a) and the 2n bits for storing the output of BWT in another array in (b): this means that during the BWT construction, we can fit almost 33 times more text using the same main memory size, thus eliminating the usage of the slower external memory for this time-consuming task in these cases.

From the combinatorial point of view, the in-place BWT is an interesting question to solve on strings. There are space saving approaches storing the suffix sorting in compressed form [START_REF] Wing-Kai Hon | A space and time efficient algorithm for constructing compressed suffix arrays[END_REF][START_REF] Joong | Alphabet-independent linear-time construction of compressed suffix arrays using o(n log n)-bit working space[END_REF][START_REF] Wing-Kai Hon | Breaking a time-and-space barrier in constructing full-text indices[END_REF][START_REF] Okanohara | A linear-time Burrows-Wheeler transform using induced sorting[END_REF] or only partially at a time [START_REF] Kärkkäinen | Fast BWT in small space by blockwise suffix sorting[END_REF], but none of them provides an in-place algorithm. In-place selection and sorting does not seem to help either [START_REF] Dobkin | Optimal time minimal space selection algorithms[END_REF][START_REF] Gianni Franceschini | In-place suffix sorting, automata, languages and programming[END_REF][START_REF] Hoare | Algorithm 65: find[END_REF][START_REF] Vinay | Selection from read-only memory and sorting with minimum data movement[END_REF][START_REF] Venkatesh Raman | Improved upper bounds for time-space trade-offs for selection[END_REF]. It is well known that in-place sorting requires the same comparison cost of (n log n) as in standard sorting. But for the BWT, we only know its comparison cost of (n log n) for the standard construction. As far as we know, no result is known for the in-place construction of BWT: a naive solution is not that simple, even if it results in exponential time. Indeed, any movement of a character T [ j] to another position inside T at least changes the content of its suffixes T i for 0 ≤ i ≤ j, making the algorithmic flavor of this problem different from that of in-place sorting n elements.

The above discussion suggests that a careful orchestration of the movement of the characters inside T is needed to avoid losing the content of some suffixes before they contribute to the BWT. Our idea is to define a sequence of transformations B 0 , B 1 , . . . , B n , where B 0 is the input text T and B n is the final BWT of T . For 1 ≤ ≤ n, we have that B is the BWT of the last characters in T and is computed from B -1 (re)using just O (1) extra memory cells. We think that this sequence of transformations could be of independent interest for the community of string algorithms, and some of the combinatorial properties that we use can be found in [START_REF] Belazzougui | Linear time construction of compressed text indices in compact space[END_REF][START_REF] Wing-Kai Hon | Breaking a time-and-space barrier in constructing full-text indices[END_REF][START_REF] Lam | Compressed indexing and local alignment of DNA[END_REF][START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler transform[END_REF].

In this paper we propose an O (n 2 )-time approach that builds the above sequence of transformations using four integer variables and one character variable, taking O (n) time per transformation in the worst case. The resulting in-place algorithm is simple and can be easily encoded in few lines of C code or similar programming languages. However we do not claim any practicality of our solution due to its quadratic worst-case cost. Our contribution is that it could lay out the path towards faster methods for the space-efficient computation of the BWT: any method to compute

B k from B k-1 in t(n) time (re)using s(n) space, would lead to a construction of the BWT in O (n • t(n)) time using O (1 + s(n)) space.
To this end, it is worth noting that the inputs for BWT are typically large and a fast algorithm that is in-place or uses very low additional memory, would be relevant in practice.

Our theoretical study has also an impact on the practical algorithm design. A natural question is what we get if we allow for some extra space. We prove that using O (k • σ k ) additional space for any given parameter k ≤ n, where σ k ≤ min{| |, k} is the maximum number of distinct characters found among k consecutive positions in T , we can compute the BWT (and

its inverse) of a text of n characters in O ((n 2 /k + n) log k)
time in the comparison model. We observe that σ k is practically a constant in many applications. The practical implications of this trade-off in space versus time can be appreciated by observing that, for any arbitrarily small > 0, we can obtain the BWT of a text of n bytes in O (n -1 log n) time using just n extra bytes for a constant-size alphabet. This is useful when the text occupies a great part of the available main memory, and only n free cells are available. This avoids using external-memory algorithms, which are clearly slower as I/O access takes several orders of magnitudes more time than main memory access. The paper is organized as follows. We describe how to perform the in-place BWT in Section 2. We then discuss how to invert the BWT, so as to obtain the original text T , in Section 3. We illustrate the trade-off between space and time in Section 4. Finally, we draw some conclusions in Section 5.

In-place BWT

Given the input text T = T [0 . . n -1] where T [n -1] = $, moving a single character inside T can change the content of many suffixes. The idea to circumvent this difficulty without using storage for the suffix sort is to proceed by induction from right to left in T , while maintaining the BWT of the current suffix T s , denoted by BWT(T s ). We assume 0 ≤ s ≤ n -3, since the last two suffixes of T are equal to their respective BWT.

To compute BWT(T s ), suppose that BWT(T s+1 ) has been already computed and stored in the last positions of T , i.e.

T [s + 1 . . n -1]. Consider the current character c = T [s]: if we look at the content of T [s . . n -1],
we no longer find T s , but the character c followed by the permutation BWT(T s+1 ) of T s+1 . Nevertheless, we still have enough information as we will show in the proof of Theorem 1 that the position of $ inside BWT(T s+1 ) is related to the rank of T s+1 among the suffixes T s+1 , . . . , T n-1 in lexicographic order. We exploit this fact in the following steps (see Fig. 1).

1. Find the position p of the $ in T [s + 1 . . n -1]: note that ps is the (local) rank of the suffix T s+1 that originally was starting at position s + 1. 2. Find the rank r of the suffix T s (originally in position s). Using character c, scan T [s + 1 . . n -1] and compute the sum of two counts: how many characters are strictly smaller than c, and how many occurrences of c appear in T [s + 1 . . p] (and add s as an offset to obtain r).

3. Store c into T [p] (thus replacing the $). 4. Insert the character $ in T [r] by shifting T [s + 1 . . r] by one position to the left, so as to occupy positions s, . . . , r -1 of T .

The C code reported in Fig. 2 implements Steps 1-4, where END_MARKER denotes $. For example, consider T = mississippi$ and s = 4, where we use capital letters to denote the BWT partially built on the last positions of T . Suppose that we have already computed the BWT for the last 7 characters in T , namely, we have missiIPSPIS$. We then have p = 11 and, since there is one character ($) smaller than c = i, and two characters that are equal to c and occur before position p, we have r = s + 3 = 7. This means that we have to replace $ by c and shift IPS by one position left so as to insert $ in position r. The next configuration is missIPS$PISI, which ends with the BWT of T s .

Theorem 1. Given a text T of n characters, we can compute its Burrows-Wheeler Transform (BWT) in O (n 2 ) time in the comparison model using O (1) additional memory cells.

Proof. We prove first the correctness. Let T be the input text and T be its modification at a generic iteration s, where

0 ≤ s ≤ n -3. Note that T [0 . . s] = T [0 . . s] while T [s + 1 . . n -1] = BWT(T [s + 1 . . n -1]
). By induction, the position p of $ in T [s + 1 . . n -1] indicates the rank ps of T s+1 among the suffixes in {T s+1 , T s+2 , . . . , T n-1 } in lexicographic order. The base case for T n-2 and T n-1 is trivially satisfied. Hence, we show how to preserve this property for 0 ≤ s ≤ n -3.

First note that the character c = T [s] goes in position p, since it precedes T s+1 inside T . Next, we have to find the new position r for T s , so that r -(s -1) is its rank among the suffixes in S = {T s , T s+1 , . . . , T n-1 } in lexicographic order. First count how many characters smaller than c occur in T [p . . n -1]: there are as many suffixes in S that are smaller than T s since their first character is smaller than c. To this quantity, add the number of occurrences of c in T [s + 1 . . p -1]: these suffixes are also smaller since they start with c but have rank smaller than p, i.e. the rank of T s+1 . In this way, we discover how many suffixes are smaller than T s in S: inserting $ in the corresponding location r of T , by shifting the characters in T [s + 1 . . r] to the left, which thus occupy positions s, . . . , r -1 (see Fig. 1), we maintain the induction. Hence,

T [0 . . s -1] = T [0 . . s -1] and T [s . . n -1] = BWT(T [s . . n -1]). When s = 0, we obtain the BWT of T .
As for the complexity, note that each of the n -2 iterations requires O (n) time, since it can be implemented by O (1) scans of T [s . .n -1]. This gives a total cost of O (n 2 ). We use four integer variables (i, p, r, s) and one character variable (c) in the C code shown in Fig. 2, and thus we need O (1) memory cells for the local variables. 2

Inverting the BWT

Reversing the permutation performed by the in-place BWT is called inverting the BWT. Initially we have the BWT of the original input text T , denoted BWT(T ). We want to invert the latter by permuting its characters. Thus we reverse the approach described in Section 2. We maintain the invariant that there is a pointer L to a certain position in the input buffer storing BWT(T ) so that, at any time, (a) the prefix of the buffer to the left of L stores the prefix of T obtained so far by the inverting process and (b) the remaining suffix of the buffer (pointed by L till the end of the input buffer) stores the portion of the BWT still to be inverted. For the sake of notation, we identify L with the entire suffix of the input buffer that still has to be inverted. Under this invariant, which is initially true by setting L to the beginning of the input buffer for BWT(T ), we proceed as follows. We find the position p of $ in L, and then select the pth character in the alphabet order in the multiset given by the characters of L. Stability is needed, since equal characters should be considered in the order of their appearance in L, as detailed below.

( i = 0; i < n; i++ ){ if (L[i] < c) count++; } / * step 3 * / f = p -count; q = -1; while ( f > 0 ){ q++; if ( L[ q ] == c ) f--; } / * step 4 * / L[ q ] = END_MARKER; for ( i = p-1; i > 0; i--){ L[i] = L[i-1]; } / * step 5 * / L[0] = c; L++; n--; / * step 1 * / if (p-1 > q)
p = q+1; / * also the new END_MARKER has been shifted * / else p = q; } } Fig. 3. Reverting the permutation of the inverse BWT.

1. Find the position p of the $ in L, and increment p (since array indexing starts from 0).

2. Let select be a selection algorithm that works on read-only input, i.e., it does not move elements around while finding the pth smallest element. Using select on L, select the pth character c in the multiset of the characters of L or, equivalently, the pth character in the sorted list of characters of L. 3. Let q denote the position inside L of the f th occurrence of c, which we hit in a stable fashion when finding c in L.

Here f is the difference between p and the number of characters c of L such that c < c. [START_REF] Burrows | A block-sorting lossless data compression algorithm[END_REF]. Replace the occurrence of c at position q by $, and remove the old occurrence of $ by shifting to the right the first p characters of L. 5. At this point, the first position in L is free: store the character c in it, and shorten L by one character at the beginning (i.e. advance the pointer L by one position towards the end of the input buffer).

The C code in Fig. 3 implements Steps 1-5 above, where END_MARKER denotes $. Note that it is a bit longer than the code for the in-place BWT in Fig. 2. As it can be seen below, the original text is reconstructed from left to right as a prefix of increasing length (indicated with small letters).

IPSSM$PISSII

→ mIPSS$PISSII → miIPSSPISSI$ → misIPSSPIS$I → missIPS$PISI → missiIPSPIS$ → missisIPSPI$ → mississIP$PI → mississiIPP$ → mississipIP$ → mississippI$ → mississippi$
The proof of correctness proceeds along the same lines as in the proof of Theorem 1, since we are reversing the procedure described there. As for the complexity, each of the n -2 iterations is dominated by the cost of select.

Let t s (n) be the time complexity in the comparison model and s s (n) be the space complexity required by select. Using the result in [START_REF] Vinay | Selection from read-only memory and sorting with minimum data movement[END_REF], we have t s (n) = O (n 1+ ) in the worst case for any fixed small constant > 0 with s s (n) = O (1), and we have t s (n) = O (n 1+ ) = O (n log log n) on the average (which meets the randomized lower bound in [START_REF] Chan | Comparison-based time-space lower bounds for selection[END_REF]), with s s (n) = O [START_REF] Adjeroh | The Burrows-Wheeler Transform: Data Compression, Suffix Arrays, and Pattern Matching[END_REF].

We can state the complexity in general terms.

Theorem 2. Let t s (n) be the time complexity in the comparison model and s s (n) be the space complexity required by select. Given the BWT of a text T of n characters, we can recover T by permuting the BWT (also known as inverse BWT) in O (n • t s (n)) time in the comparison model using O (1 + s s (n)) additional memory cells.

We give some examples of the bounds that can be attained with Theorem 2.

Corollary 1.

Given the BWT of a text T of n characters, we can recover T by inverting the BWT in O (n 2+ ) time in the worst case, or O (n 2 log log n) time on the average, in the comparison model using O [START_REF] Adjeroh | The Burrows-Wheeler Transform: Data Compression, Suffix Arrays, and Pattern Matching[END_REF] additional memory cells.

Using slightly more additional space than a constant-literally speaking, the algorithm is no more in-place-and the result in [START_REF] Venkatesh Raman | Improved upper bounds for time-space trade-offs for selection[END_REF], where t s (n) = O (n(log n) 2 ) and s s (n) = O (log n), we derive the following.

Corollary 2. Given the BWT of a text T of n characters, we can recover T by inverting the BWT in O ((n log n) 2 ) time in the comparison model using O (log n) additional memory cells.

Finally, for the special case in which the alphabet of the distinct characters in T is of constant size (as in DNA and ASCII texts), we obtain an improved bound since select can be immediately implemented by a simple scheme that employs

O (| |) = O (1) counters.

Corollary 3. Given the BWT of a text T of n characters drawn from a constant-size alphabet, we can recover T by inverting the BWT

in O (n 2 ) time in the comparison model using O (1) additional memory cells.

Practical trade-off between space and time

The inplace algorithms described so far have the drawback of requiring (n 2 ) time, which make them unfeasible for long texts. A natural question is how much the latter bound can be improved using extra space. For example, using the dynamic wavelet tree data structure [START_REF] Navarro | Optimal dynamic sequence representations[END_REF] in additional O (n + | | log n) bits of space, we can maintain the BWT through insertion and deletion operations of individual symbols, supporting rank and select operations, with a cost of O (log n/ log log n) time per operation. Using the latter data structure, our algorithms in Sections 2 and 3 would give a bound of O (n log n) time with additional O (n + | | log n) bits of space besides that needed for storing the n characters of the input text T . 3 However the resulting solution is not very practical as the data structure in [START_REF] Navarro | Optimal dynamic sequence representations[END_REF] is quite sophisticated. We show next how to smoothly adapt our algorithms in Sections 2 and 3 to a situation where extra memory is allowed, producing some trade-off solutions that are amenable for implementation with a flexible parameter k for the additional space.

Theorem 3. Given a text T of n characters, we can compute its Burrows-Wheeler Transform (BWT) and its inverse in O ((n

2 /k + n) log k) time in the comparison model using O (k • σ k ) additional space,
where σ k ≤ min{| |, k} is the maximum number of distinct characters found among k consecutive positions in T .

To appreciate the bound in Theorem 3 from a practical point of view, consider the situation in which the text T occupies a great part of the available memory, and the remaining free cells are a constant fraction of the text size. Our algorithm takes O (n log n) time by fixing k to be a suitable fraction of n. This avoids to use external-memory algorithms, which are clearly slower as I/O access takes several order of magnitudes with respect to main memory access. In general, if the available memory size is M, we obtain the following result by setting k = (Mn).

Corollary 4.

Let M be the number of available cells in main memory. Given a text T of n < M characters over a constant alphabet, we can compute its Burrows-Wheeler Transform (BWT) and its inverse in O (( n 2 M-n + n) log n) time in the comparison model using ≤ M total memory cells including those containing T . When M ≥ (1 + )n for a constant > 0, this gives O (n log n) time using just n additional cells.

The idea to prove Theorem 3 is to have n/k batches. Each batch simulates k consecutive iterations in the external for loop on s in Fig. 2, taking O ((n + k) log k) time and using O (k • σ k ) space as follows. 3 This theoretical solution has been suggested by Rossano Venturini (private communication).

Base case. Let s 1 be the largest multiple of k that is smaller than n. We can compute BWT(T [s 1 . . n -1]) and store it in

T [s 1 . . n -1] in O (k log k) time and O (k) additional space by observing that |T [s 1 . . n -1]| < k.
Inductive case. Suppose that BWT(T [s 1 . . n -1]) has been already stored in T [s 1 . . n -1], where s 1 > 0 is now a generic multiple of k. Letting s 0 = s 1k, we want to show how to store BWT(T [s 0 . .

n -1]) in T [s 0 . . n -1] in O ((n + k) log k) time using O (k • σ k ) additional space.
Since we have one base case and ≤ n/k inductive cases, the final cost will be O (k

log k + (n/k) • (n + k) log k) = O ((n 2 /k + n) log k) time using O (k • σ k ) additional space,
as stated in Theorem 3.

Inductive case

We can abstract the problem for a string X (i.e. T [s 0 . . n -1]) of length m, where the characters in X[k . . m -1] are already permuted according to their BWT, and the characters in X[0 . . k -1] are still in their original order. We want to compute BWT( X) in O ((m

+ k) log k) time using O (k • σ k ) additional space. Let Z denote X[k . . m -1]
where the $ character is virtually removed from its position, say j. Hence Z is of length mk -1 and the pair $, j is a breakpoint for Z . In general, a breakpoint is a pair c, j such that c is virtually occupying position j of Z : if two or more breakpoints claim the same position j, there should be a relative order among them.

Our goal is to compute the k + 1 breakpoints for Z so that (a) their characters are those in X[0 . . k -1] plus $, and (b) flattening Z and these breakpoints correctly produces BWT( X) as follows. Given Z and an ordered list of k + 1 break- points B = c 0 , j 0 , . . . , c k , j k , where 0 ≤ j 0 ≤ • • • ≤ j k ≤ |Z |, flattening Z and B produces a string with the characters of Z suitably shifted to the left to make room for the characters in the breakpoints of B as follows. We scan Z starting with j = 0: if j = j r for the breakpoint c r , j r at the beginning of B, we output c r and remove c r , j r from B; else ( j = j r ), we output Z [ j] and increase j. (If c r = $, we do not really output it, but we retain its position j r for the next batch.) The computation ends when Z has been completely scanned and the list B has been emptied. The required time is O (m log k) and the computation can be performed using O (k • σ k ) additional space.

For this we need the following auxiliary data structures for string Z , which require O (k • σ k ) additional space. (Note that Z and X require just O (1) space as they originate from T .)

1. Static array C of σ k ≤ k entries, where α 1 < • • • < α σ k are the distinct characters in X[0 . .k -1]: entry C [i] is the number of positions j in Z such that Z [ j] < α i . 2.
Static rank data structure R 1 on Z supporting queries that, for an integer j and a character α i , report how many positions j satisfy Z [ j] = α i and j ≤ j .

3. Dynamic list B of breakpoints, initially containing only the pair $, j .

4. Dynamic rank data structure R 2 on B supporting queries that, for a character α i , report how many breakpoints c, l to the left of $, j in B satisfy c = α i .

As it is clear, we want to populate the list B by simulating the algorithm described in Section 2. Namely, we want to find the breakpoint of character X[s ] for s = k -1, k -2, . . . , 0 in Z .

Consider the breakpoint $, j , which exists in B by construction. Let α i = X[s ], and r = r 0 + r 1 + r 2 be sum of three quantities: the number r 0 of positions j in Z such that Z [ j ] < α i , the number r 1 of positions j in Z such that Z [ j ] = α i and j ≤ j, and the number r 2 of breakpoints c, l that are to the left of $, j in B and have c = α i . Note that we can compute r 0 using entry C [i] in point 1, r 1 using the data structure R 1 in point 2, and r 2 using the data structures B and R 2 in points 3-4.

Fact 1. The value of r is the rank of

X[s . . m -1] among X[s + 1 . . m -1], X[s + 2 . . m -1], . . . , X[m -1 . . m -1] in lexicographic order. Proof. It follows from the fact that X ≡ T [s 0 . . n -1] and X[s + d . . m -1] ≡ T s+d , where s = s 0 + s and d ≥ 0. 2
After computing r , we replace the breakpoint $, j by c, j , and create a new breakpoint $, r to be inserted in B, updating the data structures in points 3-4 accordingly. Proof. We scan X[0 . . k -1] and find the distinct characters α 1 < • • • < α σ k . We then perform a scan of Z to store in C [START_REF] Adjeroh | The Burrows-Wheeler Transform: Data Compression, Suffix Arrays, and Pattern Matching[END_REF] the number of positions j such that Z [ j] < α 1 and, for i > 1, to store in C [i] the number of positions j such that α i-1 ≤ Z [ j] < α i . We then store in C its prefix sums, thus giving the wanted C . Time is O (m log σ k ) since during the scan we perform a binary search among the σ k characters. Space is O (σ k ) by definition of C . 2 Lemma 2. The static data structure R 1 can be stored in O (k • σ k ) space and built in O (m log σ k ) time, so that each query requires O (log σ k + m/k) time.

Proof. We store R 1 as a collection of σ k arrays R i 1 , for 1 ≤ i ≤ σ k . Entry R i 1 [h], for 0 ≤ h ≤ k, stores the number of occur- rences of character α i in the hth segment of m/k consecutive positions in Z : namely, R i

1 [0] = 0 and, for h ≥ 1, R i 1 [h] stores the number of positions j such that Z [ j] = α i and (m/k) • (h -1) ≤ j ≤ max{m -1, (m/k)h -1}.
After initializing the entries in all the arrays to zero, their correct value can be set with a single scan of Z in O (m log σ k ) time. 4 After that, we perform a post-processing and store in each R i 1 its prefix sums: in this way, R i 1 [h] stores the number of positions j such that Z [ j] = α i and 0 ≤ j ≤ max{m -1, (m/k)h -1}. For a query with a character c and an integer j , it takes O (log σ k ) time to establish that c = α i for a certain i, and O (1) time to find the largest h such that (m/k)h < j . The correct answer for the query is then given by the sum of the content of R i Proof. We handle the list B = c 0 , j 0 , . . . , c r , j r , where 0 ≤ r ≤ k, and the data structures R 2 together. In particular, R 2 is the wavelet tree of height O (log r) built on the sequence c 0 • • • c r of characters taken from B. The query for α i can be performed as a count query of characters α i in c 0 • • • c d , where d ≤ r and c d = $ (e.g. [START_REF] Navarro | Wavelet trees for all[END_REF]). Each update (insertion, deletion, replacement) can be handled in O (log r + log k) = O (log k) time [11,26]. 2

We now have all the ingredients to prove the time and space bounds for computing the BWT as stated in Theorem 3. By the inductive scheme and Fact 1, the computation is correctly performed. As for the time bounds, we use Lemmas 1-3. Note that to invert the BWT, we can now implement the algorithm described in Section 3 in a simpler way, since R 1 and R 2 support also the selection of the pth symbol c = α i . Moreover, flattening Z and B removes the positions j stored in the pairs in B from Z , as this simulates the shift of some characters of Z to the right. The time and space analysis is similar to that of computing the BWT.

Conclusions

We presented an in-place BWT construction taking O (n 2 ) time in the comparison model. It would be interesting to improve this bound. Note that the while loop in our in-place BWT can be avoided using O (| |) space, where is the alphabet of characters occurring in T . Time can be further reduced to O (n 2 / log | | n) by packing characters but this is still not useful for large text collections.

We do not know whether a lower bound better than (n log n) holds for the problem in the comparison model since the space is very constrained. This is an interesting question to investigate.

On the practical side, our in-place algorithms can be adapted to provide a trade-off between space and time, when O (k • σ k ) extra memory cells are allowed, providing an O ((n 2 /k + n) log k)-time algorithm to obtain (and invert) the BWT.

For example in real systems where the alphabet size is a constant, for any arbitrarily small > 0, the BWT of a text of n bytes can be computed in O (n -1 log n) time using just n extra bytes.

Fig. 1 .

 1 Fig. 1. An illustration of Steps 1-4 of the in-place construction of the BWT.

Fig. 2 .

 2 Fig. 2. In-place construction of BWT.

  void inplaceIBWT( unsigned char L[ ], int n ){ int f, i, p, q, count;

Lemma 1 .

 1 The static array C can be stored in O (k • σ k ) space and built in O (m log k) time.

1 2 Lemma 3 .

 23 [h] and the number of α i 's found in Z [(m/k)h . . j ] by its direct inspection.Scanning the latter takes O (m/k) time by definition of R i 1 . The dynamic list B and the dynamic data structure R 2 can be stored in O (k) space and built in O (k log k) time, so that each query and each update requires O (log k) time.

Table 1 BWT

 1 L for the text T = mississippi$ and its relation with suffix sort.

	Cyclic shifts	Sorted cyclic shifts		Suffixes
			L	i	T i
	mississippi$	$ mississipp	i	11	$
	$mississippi	i$mississip	p	10	i$
	i$mississipp	ippi$missis	s	7	ippi$
	pi$mississip	issippi$mis	s	4	issippi$
	ppi$mississi	ississippi$	m	1	ississippi$
	ippi$mississ	mississippi	$	0	mississippi$
	sippi$missis	pi$mississi	p	9	pi$
	ssippi$missi	ppi$mississ	i	8	ppi$
	issippi$miss	sippi$missi	s	6	sippi$
	sissippi$mis	sissippi$mi	s	3	sissippi$
	ssissippi$mi	ssippi$miss	i	5	ssippi$
	ississippi$m	ssissippi$m	i	2	ssissippi$

As is standard in many string algorithms, we assume that any two characters in can only be compared and this takes O (1) time. Hence, comparing characterwise any two suffixes may require O (n) time in the worst case.

In C code, we would declare T as unsigned char T[n] and use this storage plus O (1) local variables of constant size. A more formal model would say that each memory cells hosts a character from and so an integer of log n bits requires log | | n cells. We prefer to keep it simpler and say that an extra cell can contain an integer.

If k • σ k = w(n log σ k ), it is a standard trick to allocate O (k • σ k ) memory and initialize only what is needed in O (m log σ k ) time[START_REF] Aho | The Design and Analysis of Computer Algorithms[END_REF].

Acknowledgements

The second author is grateful to Gianni Franceschini for some preliminary discussions and to Venkatesh Raman for pointing out the results in [23,28] and Rossano Venturini for his comments and indicating the result in [3]. We also thank the anonymous reviewers for their careful reading of our manuscript.


Simone Hemm

François Caire

Jerome Coste

François Vassal

Christophe Nuti

Philippe Derost

Lemlih Ouchchane

Laurent Sarry

Franck Durif

Jean-Jacques Lemaire
email: jjlemaire@chu-clermontferrand.fr

Jérôme Coste

Ouchchane Lemlih

Postoperative control in deep brain stimulation of the subthalamic region: the contact membership concept

Keywords: STN-DBS, Direct targeting, MRI, Postoperative control, Contact mapping

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Deep brain stimulation (DBS) of the subthalamic nucleus (STN) has proved to be a highly effective treatment of symptoms related to severe idiopathic Parkinson's disease (PD) [START_REF] Hamani | Bilateral subthalamic nucleus stimulation for Parkinson's disease: a systematic review of the clinical literature[END_REF][START_REF] Krack | Five-year follow-up of bilateral stimulation of the subthalamic nucleus in advanced Parkinson's disease[END_REF][START_REF] Krause | Chronic inhibition of the subthalamic nucleus in Parkinson's disease[END_REF][START_REF] Rodriguez-Oroz | Efficacy of deep brain stimulation of the subthalamic nucleus in Parkinson's disease 4years after surgery: double blind and open label evaluation[END_REF].

However, the anatomic structures affected by stimulation, and therefore the underlying mechanisms of action of DBS, have not yet been determined [START_REF] Benazzouz | High-frequency stimulation of both zona incerta and subthalamic nucleus induces a similar normalization of basal ganglia metabolic activity in experimental parkinsonism[END_REF][START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Henderson | Postmortem analysis of bilateral subthalamic electrode implants in Parkinson's disease[END_REF][START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF][START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Zincone | Physiologic study of the subthalamic volume[END_REF]. The success of DBS electrode implantation resides essentially in the relief of symptoms, but the position of electrode contacts needs to be determined in terms of anatomic position. More and more groups advocate postoperative imaging (magnetic resonance imaging MRI/computer tomography CT) to check for absence of hemorrhage and determine the final electrode position in comparison with the planned one [START_REF] Duffner | Relevance of image fusion for target point determination in functional neurosurgery[END_REF][START_REF] Ferroli | A simple method to assess accuracy of deep brain stimulation electrode placement: pre-operative stereotactic CT postoperative MR image fusion[END_REF][START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF]. Only a few of these groups determine the contact positions and correlate them with anatomy [START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF][START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Nowinski | A probabilistic functional atlas of the human subthalamic nucleus[END_REF][START_REF] Starr | Implantation of deep brain stimulators into the subthalamic nucleus: technical approach and magnetic resonance imaging-verified lead locations[END_REF][START_REF] Vayssière | Deep brain stimulation for dystonia confirms a somatotopic organization in the Globus Pallidus Internus[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF]. Still fewer groups base this analysis on the patient's MRI alone without atlas projection [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Vayssière | Deep brain stimulation for dystonia confirms a somatotopic organization in the Globus Pallidus Internus[END_REF]. Nearly all of the groups performing a detailed anatomic study localize the anatomic position of centers of the electrode contacts. They do not consider the contact in its dimension [START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Nowinski | A probabilistic functional atlas of the human subthalamic nucleus[END_REF][START_REF] Starr | Implantation of deep brain stimulators into the subthalamic nucleus: technical approach and magnetic resonance imaging-verified lead locations[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF][START_REF] Zonenshayn | Location of the active contact within the subthalamic nucleus (STN) in the treatment of idiopathic Parkinson's disease[END_REF].

The work reported here describes the postoperative position analysis of the clinically effective contacts, i.e., contacts used for chronic stimulation, based on individual patient anatomy. In accordance with our targeting concept (preoperative direct targeting linked to manual structure segmentation), it is based on the visual analysis of pre-and postoperative MR images. As the electrode size (diameter = 1.27mm; height = 1.5mm) is significant compared with the functional targets, it is often difficult to assign a given contact to a single brain structure. We therefore propose a new semi-quantitative approach, the contact membership concept, which takes into account the extent of involvement for each anatomic structure, i.e., the initial target, the STN, and its neighborhood. The classical binary approach, i.e., where the center of the contact is [START_REF] Hamani | Bilateral subthalamic nucleus stimulation for Parkinson's disease: a systematic review of the clinical literature[END_REF] or is not (0) in a given structure, differs from our membership concept, where the whole contact is related to different structures: anatomic membership of a contact; ordinal values from 0 to 1, weighting the geometric relationships with the structures.

MATERIALS AND METHODS

Patients

Ten patients suffering from severe idiopathic Parkinson's disease and treated by bilateral monopolar STN-high frequency stimulation according to the guidelines of the French National Health Agency for Accreditation and Evaluation were included in this study [5 males, 5 females; age 61.8 ± 8.1years, duration of disease 12.7 ± 5.2years (mean ± SD)]. Their clinical improvement was evaluated using the unified Parkinson's disease rating Scale (UPDRS) [START_REF] Fahn | Unified Parkinson's disease rating scale[END_REF]. It was determined in OFF-drug condition (for at least 12h) with a mean voltage of 2.6 ± 0.5V (frequency 136.0 ± 12.3Hz, pulse width 63.0 ± 9.2 μsec). Global improvement after 3 months of follow-up was 50.9 ± 12.0% (range 33.7-71.4%). The lateralized UPDRS III mean motor subscore improvement (sum of items 20-26 for each body side, resulting in a score range of 0-36), was 61.8 ± 16.1% (range 40.0-100.0%). This measure reflects the degree of tremor, rigidity, and bradykinesia in the left and right extremities, excluding axial symptoms, and has already been proposed by other authors for the evaluation of the contralateral effects of left and right electrodes [START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF]. Contacts resulting in the best alleviation in rigidity with the lowest voltage and without significant side effects were considered as optimal and selected for chronic stimulation as the effective contacts. The total number of effective contacts available for analysis was 20

(pooling right and left sides).

Surgical procedure

Stereotactic electrode implantation (Leksell G frame, Elekta, Sweden) was performed under local anesthesia in a 2 day procedure. A preoperative T2-weighted stereotactic MRI was obtained and the STN and neighboring structures were identified and segmented manually using stereotactic software (iPlan, Brainlab, Germany). On the first day we performed a White Matter Attenuated Inversion Recovery (WAIR) T2-weighted sequence (Table 1) to visualize both the stereotactic markers and the subthalamic anatomy of nuclei and bundles [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Magnotta | Visualization of subthalamic nuclei with cortex attenuated inversion recovery MR imaging[END_REF]. Several structures of the subthalamic region were outlined manually on the coronal plane (see example Fig. 3): STN, zona incerta (ZI), substantia nigra (SN), red nucleus (RN), nucleus ansa lenticularis (AL), the H1 (thalamic fascicule) or H2 (lenticular fascicule) Forel's fields (FF), the mamillo-thalamic bundle (MTB), the ventral thalamus (Thal) and the periventricular region (PVR). The recognition of the anatomy was achieved through a step-by-step approach starting from the easiest identification of simple shaped and located structures such as SN, RN and STN [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF]. According to the relative position and with the help of anatomic and stereotactic books and in-house 4.7-tesla (T) MRI anatomy software an accurate 3D multiplane analysis was then possible [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Derost | Is subthalamic nucleus deep brain stimulation (DBS-STN) appropriate to manage severe Parkinson disease in an elderly population[END_REF][START_REF] Lemaire | Contribution of 4.7 tesla MRI in the analysis of the MRI anatomy of the human subthalamic area[END_REF].

Table 1. Technical characteristics of MRI sequences using a head coil: pre-operative, stereotactic (stereotactic frame in place), T2-weighted; post-operative, non-stereotactic, T1weighted, Specific absorption rate = 0.1W/kg.

On the second day, after repositioning of the stereotactic frame (Repositioning kit, Elekta, Sweden), the optimal target was defined under local anesthesia based on the anatomic information, intraoperative microelectrode recording (MER) and the clinical benefit and adverse effects assessment following a gradual acute stimulation. A quadripolar electrode (DBS 3389, Medtronic, USA) was placed during the same procedure with the center of contact 1 (right electrode) or 5 (left electrode), located just above the distal contact (0, right electrode or 4, left electrode) on this optimal point. Anysignificant brain shift and electrode misplacement (checked by intra and postoperative radiologic controls) were avoided during the bilateral electrode implantation thanks to precautions: (1) the patient operated in a recumbent position in order to maintain a positive intracranial pressure, (2) the "closed dura matter" approach allowed by water-air tight guides (Bengun, Immi, France), (3) double obliquity trajectories in order to avoid ventricles, (4) two exploration tracts (two long guide-tubes) and ( 5) brain maintained by one guide-tube (on the trajectory which was not selected after the intraoperative tests) while the DBS electrode is placed coaxially along the selected tract. To check the final electrode position, a stereotactic X-ray control (based on Leksell Coordinates Indicator Box) was performed directly after electrode placement, and a T1-weighted, non-stereotactic 3D-MPRAGE MRI control in the axial plane (Table 1) within 2days. A neuropacemaker (Kinetra, Medtronic, USA) was implanted in the abdominal region under general anesthesia some days later.

Contact localization

Following our surgical procedure, the final electrode position could be determined via the postoperative MRI control. To collect preoperative stereotactic data, the segmented structures and the planned trajectories on the postoperative MRI, we merged postoperative (non-stereotactic) MRI with preoperative (stereotactic) MRI using the mutual information algorithm (iPlan, BrainLAB, Germany) [START_REF] Maintz | A survey of medical image regist ration[END_REF][START_REF] Wells | Multi-modal volume registration by maximization of mutual information[END_REF]. The final electrode position was determined based on the black electrode artifact on the postoperative MRI. As shown elsewhere, the real electrode position can be assigned to the middle of the artifact [START_REF] Pollo | Magnetic resonance artifact induced by the electrode Activa 3389: an in vitro and in vivo study[END_REF]. The position of each contact was calculated starting from the middle of the artifact (the center position enclosed by the four contacts of the electrode) and based on the known theoretical electrode dimensions (DBS 3389: 7.5mm between the ends of the proximal and distal contacts, 0.5mm between contacts). The anatomic relationships of the contacts were determined on reconstructed images along the electrode to optimize the interpretation of anatomic position.

Definition of the anatomic position

The analysis of the anatomic position of a given contact was performed using three views along the electrode (called "pseudo" coronal, sagittal and axial). We considered the contact position in relationship with the neighboring structures, i.e., the STN, ZI, SN, AL, Forel's fields H1 and H2, and the ventral thalamus. The STN was subdivided into three different parts, a lateral superior part (lat-sup STN), an intermediate part (int STN) and a medial inferior part (med-inf STN).

The identification was performed in two different ways:

1. Classical approach: contact center-one structure (binary: 0/1). Only one anatomic structure was assigned to the effective contact according to the projection of the center of this contact.

2. Membership concept: contact in its size-several structures (ordinal values: from 0 to 1).

To take into account the imprecision, especially the one linked to the MRI resolution, we chose a mathematical model based on the variable "degree of membership" used in fuzzy logic introduced by Zadeh [START_REF] Zadeh | Fuzzy sets[END_REF][START_REF] Zadeh | Fuzzy logic and approximate reasoning (in memory of Gregore Moisil)[END_REF]. We subjectively assigned a weight to each anatomic structure covered by the effective contact. The weight resulted from the ordinal coding accounting for the membership degree and from the different anatomic configurations. When the contact could be clearly assigned to only one structure, the weighting was obviously set to 1. When several structures were involved, we distinguished between main and secondary structures according to the size of the contact part lying within the considered structure. The larger the number of structures involved, the lower the weighted values assigned to the structures (Table 2, Fig. 1). If there is only one structure that can be defined as a main structure but if the "contact" is shared with other structures involved in a lesser extent (see secondary structure), the weight for that main structure is set to 0.8. If there are several structures that can be defined as main structures, the weights of these structures are set to 0.5each, if they are the only structures involved. In the case of several main structures, their weights are set to 0.4 each if other (secondary) structures are also involved. A structure is referred to as secondary structure if the extent of its contact is considered as "minor". Whatever the configuration and the number of structures involved, the weight for a secondary structure is set to 0.25.

Table 2. Membership degree (from 0 to 1) of given anatomic structures according to the extent of their involvement for a considered contact; main: principal structure covered by a large part of the contact; secondary: secondary structure covered by a small part of the contact.

Data analysis

For each active contact (n=20), we classified the anatomic relationships as described above. For each structure we calculated the sum of all contact memberships, once with the binary-center approach and once with the contact membership concept. Then the percentage of involvement of each structure was calculated. Mean and standard deviations of the weightings for each structure were additionally determined. We confronted, for each patient separately, the results obtained with the membership concept with those of the classical approach based on the center of the contact.

RESULTS

The image series of the ten patients could all be analyzed for the right and left electrodes.

Accordingly, we explored the anatomic position of 20 active contacts.

With both approaches (center of contact and whole contact), the lateral part of the STN, the Forel field H1 and the zona incerta can mostly be found in the vicinity of the contact (Fig. 2). One essential difference between the approach using only the contact center and the one using the whole electrode contact is that the intermediate STN appears rarely and the AL not at all in the first approach. In the classical approach, based on the contact center, the structures most often involved were the lateral STN and the zona incerta (28.6%), followed by the Forel field H1 (19.1%). With the membership approach the frequency of appearance changed, the ZI became more frequent than the lateral STN (ZI: 32.0%, STN: 24.1%). The intermediate STN and the Forel Field H2 became almost as frequent as the H1 (H1: 13.37%, H2: 11.7%, int STN: 10.98%).

The analysis of each patient separately (Table 3, Fig. 3) confirmed that lateral STN and ZI were the structures most often involved (10 out of 20 cases, for each structure) and that in general at least one of the two structures could be found close to the contact. Forel's Fields H1 (7/20) and H2 (6/20) also appeared regularly.

The contact center approach assigns one structure to each contact (Table 3, grey cells). With the contact membership approach, we obtained such a result, i.e., only one involved structure, in only 7 out of 20 cases (Table 3, structures set to "1"). In the other 13 cases, one or two other additional structures appeared. In addition, we identified, for each contact, the structure(s) presenting the highest membership degree: in 12 cases there was only one structure with the highest membership degree (Table 3, structures set to "*"), while in 8 cases there were several principal structures with the same extent of involvement.

Table 3. Analysis of the position of effective contacts (n=20) for the ten patients (right R and left L electrodes) by the membership concept (considering contact size). The structures involved according to the binary classical approach (only considering contact center) are given in italics AL nucleus ansa lenticularis, lat-sup lateral superior part of the STN, int intermediate part of the STN a: Clinical motor improvement expressed in percentage based on the UPDRS motor score:

(preoperative scorepostoperative score)/ postoperative score × 100. b: Only one main structure has been identified by the membership concept and that this structure corresponds to the one involved according to the binary concept. Cases where there are several main structures are not included (8 out of 20 cases) Figure

Example for one patient: (Top) original coronal raw images with the outlines of anatomic structures (top row; left , overview; right, zooming, labelled structures); (Bottom) reconstructed MRI plane, pseudocoronal, along the electrode trajectory (bottom row); analysis according to the contact center ( left ) or the whole contact ( right )

. Following the first approach only one structure is involved against three in the second one. H1 H1 Forel's field (thalamic fascicle), H2 H2 Forel's field (lenticular fascicle), al ansa lenticularis complex, stn subthalamic nucleus, sn sustantia nigra, zi zona incerta.

DISCUSSION

Post-placement imaging is mainly used to check for absence of hemorrhage and also increasingly to analyze the final electrode position in DBS. The anatomic position of a given contact is less frequently sought. In most teams the anatomic position is determined on the postoperative MR images [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Vayssière | Deep brain stimulation for dystonia confirms a somatotopic organization in the Globus Pallidus Internus[END_REF][START_REF] Ondo | The North American survey of placement and adjustment strategies for deep brain stimulation[END_REF] or on X-ray controls [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF]; it is assumed that there is no significant brain shift which might cause misinterpretations [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF][START_REF] Miyagi | Brain shift: an error factor during implantation of deep brain stimulation electrodes[END_REF]. On the MR images the center of the artifact is generally considered as the middle of the four contact location [START_REF] Pollo | Magnetic resonance artifact induced by the electrode Activa 3389: an in vitro and in vivo study[END_REF] and the different electrode contacts are deduced manually [START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Starr | Implantation of deep brain stimulators into the subthalamic nucleus: technical approach and magnetic resonance imaging-verified lead locations[END_REF][START_REF] Vayssière | Deep brain stimulation for dystonia confirms a somatotopic organization in the Globus Pallidus Internus[END_REF] or automatically [START_REF] Butson | Patient-specific analysis of the volume of tissue activated during deep brain stimulation[END_REF]. Pre-and postoperative MR image matching based on the mutual information algorithm [START_REF] Wells | Multi-modal volume registration by maximization of mutual information[END_REF] has been shown to be reliable [START_REF] Maintz | A survey of medical image regist ration[END_REF] and thus can be used to analyze the contact position relative to MRI anatomy. On stereotactic X-ray images, contacts are visible separately and their positions can be determined directly, in general manually, in the stereotactic reference system [START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF][START_REF] Hamel | MRI-and skull X-ray-based approaches to evaluate the position of deep brain stimulation electrode contacts-a technical note[END_REF].

Unfortunately, it is difficult to compare results obtained with the different methods applied to localize a contact and the surrounding structures. The analysis has three different characteristics:

(1) how the contact is represented, as a point or with its real dimensions; (2) the anatomic reference used (atlas, MRI, electrophysiologic study); and (3) how the contact membership is interpreted (assignment to one or several structures, weighted or non-weighted). In the following section we discuss these three features.

Contact representation

Most groups analyze the active contact position by projecting the center of the contact on the anatomy [START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF][START_REF] Nowinski | A probabilistic functional atlas of the human subthalamic nucleus[END_REF][START_REF] Starr | Implantation of deep brain stimulators into the subthalamic nucleus: technical approach and magnetic resonance imaging-verified lead locations[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF][START_REF] Zonenshayn | Location of the active contact within the subthalamic nucleus (STN) in the treatment of idiopathic Parkinson's disease[END_REF]. Herzog et al. [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF] took into account the height of the contact (1.5mm).

None of the groups considered the real contact size. This may have been for technical reasons, such as the chosen anatomic reference (see below): when a classical stereotactic atlas is used, the contact cannot be visualized along its axis since only several slices are available, with no possibility of reorientation and interpolation as there is with MR images.

Anatomic reference

There are several possibilities to analyze the anatomic position of a contact. The first one, the stereotactic atlas, is still used by many groups [START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Nowinski | A probabilistic functional atlas of the human subthalamic nucleus[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF][START_REF] Hamel | MRI-and skull X-ray-based approaches to evaluate the position of deep brain stimulation electrode contacts-a technical note[END_REF]. The contact, represented by its center, can be projected, after normalization and with reference to the AC-PC midpoint, on the stereotactic Schaltenbrand and Bailey atlas discrete slices [START_REF] Schaltenbrand | Introduction to stereotaxis with an atlas of the human brain[END_REF]. The atlas can also be projected, after deformation, on the patient MRI.

The geometric distortion in current clinical MRI (1.5-T machines) is now approaching tenths of a millimeter [START_REF] Wang | Geometric distortion in clinical MRI systems Part II: correction using a 3D phantom[END_REF][START_REF] Wang | Geometric distortion in clinical MRI systems Part I: evaluation using a 3D phantom[END_REF]. In view of these advances in technology, another possibility for analyzing the contact position is the projection of the electrode contact directly on the patient's MRI , which is generally chosen by those groups already working at least partially with direct MRI targeting [START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF]37]. In this case, it is possible to visualize not only the center of the contact but also its height and (or) diameter.

In addition, these two methods can be combined with the results of microelectrode recording as proposed by Herzog et al. who defined the dorsolateral margin of the STN by the analysis of intra-operative microrecordings [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF].

Interpretation

Different ways to interpret the anatomic contact position have been used. Groups defining the contact by its center have often based their analysis on the classical approach, assigning one structure to the active contact [START_REF] Starr | Implantation of deep brain stimulators into the subthalamic nucleus: technical approach and magnetic resonance imaging-verified lead locations[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF] or describing the position relative to certain structures or borders [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF][START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Zonenshayn | Location of the active contact within the subthalamic nucleus (STN) in the treatment of idiopathic Parkinson's disease[END_REF]. Herzog analyzed the position of the active contact in its length relative to the distance from the electrophysiologically determined dorsal STN border [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF]. Following this approach they classified the active contacts into three groups: contacts lying above the STN (+1.5mm or above), contacts at the border zone of the STN (-1.5mm to +1.5mm) and contacts inside the STN (-1.5mm or below). Voges et al. [START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF] considered three different anatomic positions of the contact center: effective contacts lying (1) within the STN, (2) outside the STN and within subthalamic fiber tracts and (3) at the border zone between the caudal third of the STN and the subthalamic white matter (STN-fiber interface). They also established the ratio of clinical improvement and stimulation energy and they accordingly took into account the fact that the therapeutic target could be some millimeters away from the effective contact. Saint-Cyr et al. [START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF] projected their results on Schaltenbrand & Bailey atlas slices and they integrated in the analysis an average current spread of 3mm radius around the cathode. In our group, we have already proposed an MRI-based quantitative approach based on the contact center and on fuzzy logic, taking into account the uncertainties in the identification of the contact [START_REF] Caire | Subthalamic deep brain stimulation for severe idiopathic Parkinson's disease. Location study of the effective contacts[END_REF]: when a contact could be identified exactly within one structure, the value 1 was assigned; when a structure was clearly not involved, the value 0 and when there was a doubt or when several structures were involved, we assigned 0.5. In the present study, we aimed to take this approach further. The first step was to introduce the real contact size which, in our opinion, is not negligible (diameter = 1.27mm; height = 1.5mm). Our stereotactic software (iPlan, Brainlab, Munich, Germany) enabled us to plot the electrode contact in the direction of the trajectory. Following our targeting concept, we used as an anatomic basis the preoperative MR images (WAIR sequence) with the manually outlined structures identified with the help of a 4.7T 3D MRI atlas (pixel size 0.25mm 3 ) [START_REF] Lemaire | Contribution of 4.7 tesla MRI in the analysis of the MRI anatomy of the human subthalamic area[END_REF] and the advanced knowledge of deep brain structures and fiber anatomy [START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Lemaire | Brain mapping in stereotactic surgery: a brief overview from the probabilistic targeting to the patient-based anatomic mapping[END_REF]. Taking both concepts together, we were able to interpret the electrode's anatomic position in 3D (on the reconstructed planes along the electrode) without stereotactic atlas-based normalization. In this way, we directly considered the patient's anatomy without reference to anatomic landmarks and avoided the reported risks, when using atlases, of geometric imprecision due to the interndividual variability of brain shape, structure size and structure location [START_REF] Kotter | Mapping brains without coordinates[END_REF][START_REF] Vayssiere | Comparison of atlas-and magnetic resonance imaging-based stereotactic targeting of the globus pallidus internus in the performance of deep brain stimulation for treatment of dystonia[END_REF]. As in most cases the contact in its dimensions cannot be assigned to a single structure, we decided to introduce the membership concept. This approach takes into account, first, the imprecision in the determination of the contact position within the artifact, and second, the part of the contact covered by a given structure. In this way we were able to consider the extent of involvement of each structure.

Our results show that the lateral STN, the Forel's field H1 and especially the ZI are the structures most often found close to the contact, consistent with the results of other groups. In addition, our study indicates that the Forel's field H2 and the intermediate STN are certainly involved also.

Hamel reported [START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Hamel | MRI-and skull X-ray-based approaches to evaluate the position of deep brain stimulation electrode contacts-a technical note[END_REF] that most contacts were located close to the dorsal border of STN or further dorsal within the subthalamic region. The other contacts were detected within the dorsal STN. The average position of all active contacts in their population coincided with the interface formed by the antero-dorso-lateral margin of the subthalamic nucleus, the ZI, and Forel's fields H2. Some of the studies using MRI-based techniques to identify the position of the active contact reported similar results [START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF][START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF]. In the opinion of Herzog et al. [START_REF] Herzog | Most effective stimulation site in subthalamic deep brain stimulation for Parkinson's disease[END_REF], it is crucial, to achieve optimal clinical results, to position the electrode at the dorsolateral border, without locating it exclusively in the subthalamic fiber tract. This seems to be confirmed by our results and by those of other groups [START_REF] Hamel | Deep brain stimulation of the subthalamic nucleus in Parkinson's disease: evaluation of active electrode contacts[END_REF][START_REF] Lemaire | MRI anatomical mapping and direct stereotactic targeting in the Subthalamic region: relationsships between biological and anatomical data in 15 Parkinsonians[END_REF][START_REF] Zonenshayn | Location of the active contact within the subthalamic nucleus (STN) in the treatment of idiopathic Parkinson's disease[END_REF][START_REF] Lanotte | Deep brain stimulation of the subthalamic nucleus: anatomical, neurophysiological, and outcome correlations with the effects of stimulation[END_REF].

Considering the ratio of clinical improvement and the energy used, Voges et al. [START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF] reported a trend toward better results for active poles projecting onto subthalamic fiber tracts than with those placed either inside the STN or at the STN-fiber interfaces. Our comparison between the classical approach, considering only one structure per contact, and the membership concept confirms that in most cases more than one structure is present close to the contact (Table 3). The extentof involvement of all the structures changes with the membership concept. The intermediate STN and the AL rarely appear or do not appear at all with the first approach (Fig. 2).

The influence of these areas in the DBS effect is unclear, as it is not usually documented in the literature, at least not specifically. Thus we draw no conclusion concerning the ansa lenticularis involvement because of its extremely low incidence in our study. Conversely the involvement of the intermediate STN seems robust and fits with the functional segmentation of STN [START_REF] Rodriguez-Oroz | Efficacy of deep brain stimulation of the subthalamic nucleus in Parkinson's disease 4years after surgery: double blind and open label evaluation[END_REF][START_REF] Parent | Deep brain stimulation in movement disorders: stereotactic coregistration of two-dimensional electrical field modeling and magnetic resonance imaging[END_REF].

A detailed analysis based on the new concept is in progress on a larger patient population in order to correlate the anatomic contact position with the clinical result for all the contacts, not only for the therapeutic one. In the long term, we intend to integrate the stimulated volume (effective electric current) in our approach, as we agree with Butson et al. [START_REF] Butson | Patient-specific analysis of the volume of tissue activated during deep brain stimulation[END_REF], Saint-Cyr et al. [START_REF] Saint-Cyr | Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging[END_REF] and Voges et al. [START_REF] Voges | Bilateral highfrequency stimulation in the subthalamic nucleus for the treatment of Parkinson disease: correlation of therapeutic effect with anatomical electrode position[END_REF] that stimulation parameters influence the extent of stimulation and that the simple analysis of the contact position is only a first step in understanding the mechanism of action of DBS [43].

CONCLUSIONS

Postoperative imaging techniques are more and more part of routinely performed surgical protocols. We have shown that in addition to its utility as a quality control procedure, it can be used to evaluate the anatomic contact position on a routine basis. Given the size of the DBS contacts, the correlation between the clinical results and the contact position will probably be more precise if the contact is considered in its extent than if only the center of the contact is used. Also, the use of MRI allows reorientation along the contact direction and the 3D position analysis with regard to the individual anatomy. The contact membership concept allows a semiquantitative analysis, assigning a given contact to several structures with a weighted membership degree. This concept offers a promising approach for further studies of the mechanism of action of DBS in Parkinson's disease and DBS in general. Our preliminary results confirm that DBS probably acts through the three main structures, the lateral STN, ZI and FF.

Figure 1 .

 1 Figure 1. Presentation of the contact membership concept on identified structures after manual outlining on original MRI slices. a The entire region of interest around the STN . b The contact can be clearly identified within one single structure (STN:1). c The contact is largely covered by one structure (STN:0.8) and marginally by one or more other structures (H2:0.25). d Two mainstructures (STN:0.5, H2:0.5). e Main (STN:0.4,H2:0.4) and secondary (ZI:0.25) structures. f Only secondary structures (ZI:0.25,STN:0.25, H1:0.25). STN subthalamic nucleus, SN substantia nigra, RN red nucleus, H1 H1 Forel's field (thalamic fascicle), H2 H2 Forel's field (lenticular fascicle), ZI zona incerta, PRR prerubral region

Figure 2 .

 2 Figure 2. Percentage involvement of each structure according to the classical binary classification (0 or 1) based on only the contact center, and according to the contact membership concept (between 0 and 1) taking into account the whole contact size.