Stéphane Kremer

Caroline Nicolas-Ong

Thérèse Schunck

Malika Schenck

Olivier Collange

Izzie Véronique Mutschler

Jacques Namer

Jean-Louis Dietemann

Chantal Delon-Martin

Francis Schneider

Interet De

Usefulness of functional MRI associated with PET scan and evoked potentials in the evaluation of brain functions after severe brain

Keywords: brain injury, coma, MRI, functional imaging, consciousness, critical care patients

à 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 :

The prevalence, complexity and multiple medical and socio-economic consequences of severe alterations of brain functions after severe central nervous system injuries, make them a major health challenge. With few exceptions, the pathophysiology and the molecular basis of both acute brain injuries and recovery have yet to be defined, in order to promote proper and efficient treatments. Before that step, the clinical challenge will be to best evaluate injuries because both clinical examination and available tools lack accuracy not only for the diagnosis of damage but also for its prognosis. This is especially true for patients with coma of intermediate severity.

Minimally conscious state (MCS) and vegetative state (VS) are clinically defined as chronic disorders of consciousness. In the vegetative state, patients demonstrate preserved hypothalamic and autonomic functions but are unaware of themselves and of their environment, and are incapable of interaction with others, whereas in the minimally conscious state, patients have limited but discernible evidence of awareness of themselves and of their environment [START_REF] Bernat | Chronic disorders of consciousness[END_REF].

Different clinical evaluation scales of disorders of consciousness have been developed and help nosological distinction between these entities [START_REF] Giacino | The vegetative and minimally conscious states: consensus-based criteria for establishing diagnosis and prognosis[END_REF][START_REF] Giacino | The minimally conscious state: definition and diagnostic criteria[END_REF][START_REF] Majerus | Behavioral evaluation of consciousness in severe brain damage[END_REF][START_REF] Teasdale | Assessment of coma and impaired consciousness. A practical scale[END_REF][START_REF] Wijdicks | Validation of a new coma scale: The FOUR score[END_REF]. Nevertheless, it remains difficult to evaluate clinically brain functions in many patients [START_REF] Bernat | Can functional MRI detect awareness when a neurological examination does not?[END_REF]. If available such information could be of particular clinical interest, as it has been hypothesized that the presence of preserved brain functions could have a prognostic value [START_REF] Hildebrandt | Brain perfusion and VEP reactivity in occipital and parietal areas are associated to recovery from hypoxic vegetative state[END_REF][START_REF] Naccache | Auditory mismatch negativity is a good predictor of awakening in comatose patients: a fast and reliable procedure[END_REF]. The knowledge of preserved brain functions could lead to more efficient rehabilitation planning or to a better adaptation of stimulation treatment [START_REF] Giacino | Functional neuroimaging applications for assessment and rehabilitation planning in patients with disorders of consciousness[END_REF][START_REF] Laureys | How should functional imaging of patients with disorders of consciousness contribute to their clinical rehabilitation needs?[END_REF].

Electrophysiological techniques, especially evoked potentials (EP) (visual, auditory and somatosensory) are widely used in this indication [START_REF] Lew | Prognostic value of evoked and event-related potentials in moderate to severe brain injury[END_REF]. Indeed, EP are used in severely brain injured patients to predict clinical outcome, as the absence of residual brain potentials could identify those patients unlikely to recover [START_REF] Hildebrandt | Brain perfusion and VEP reactivity in occipital and parietal areas are associated to recovery from hypoxic vegetative state[END_REF][START_REF] Rothstein | The role of evoked potentials in anoxic-ischemic coma and severe brain trauma[END_REF]. The advantage of the technique is that it is easy to perform at bedside, but its major drawbacks are the frequent presence of environment -related artefacts and a gross localisation information.

Positron Emission Tomography (PET) allows evaluation of cerebral glucose metabolism whereas Single Photon Emission Computed Tomography (SPECT) allows evaluation of cerebral perfusion : They can demonstrate a decrease of global cerebral metabolism and perfusion in persistent vegetative patients. This decrease is more marked in arterial border zones in post-anoxic patients [START_REF] Beuthien-Baumann | Persistent vegetative state: evaluation of brain metabolism and brain perfusion with PET and SPECT[END_REF][START_REF] Volder | Brain glucose metabolism in postanoxic syndrome due to cardiac arrest[END_REF][START_REF] Laureys | Cerebral metabolism during vegetative state and after recovery to consciousness[END_REF][START_REF] Rudolf | Cerebral glucose metabolism in acute and persistent vegetative state[END_REF][START_REF] Schiff | Residual cerebral activity and behavioural fragments can remain in the persistently vegetative brain[END_REF], and there is a correlation between areas of relatively preserved brain metabolism and residual brain function [START_REF] Schiff | Residual cerebral activity and behavioural fragments can remain in the persistently vegetative brain[END_REF]. However, these techniques are not available at bedside everywhere and require either patient transportation or that specialized centers be located close to critical care departments.

More recently functional neuroimaging techniques (H 2 15 0 PET and functional MRI (fMRI)) have been used in clinical research studies to evaluate brain functions in patients presenting chronic disorders of consciousness [START_REF] Bernat | Chronic disorders of consciousness[END_REF][START_REF] Laureys | Functional neuroimaging in the vegetative state[END_REF][START_REF] Laureys | Brain function in coma, vegetative state, and related disorders[END_REF][START_REF] Schiff | Multimodal neuroimaging approaches to disorders of consciousness[END_REF].

Performed late after the initial injury in comatose patients, these studies have shown that auditory, nociceptive and visual passive stimulations lead to activations of primary cortices but also in some cases of secondary associative cortex. The presence of activation in higher order association areas suggests the persistence of integrative processes and may herald better prognosis [START_REF] Giacino | Functional neuroimaging applications for assessment and rehabilitation planning in patients with disorders of consciousness[END_REF][START_REF] Laureys | Brain function in coma, vegetative state, and related disorders[END_REF][START_REF] Boly | Auditory processing in severely brain injured patients: differences between the minimally conscious state and the persistent vegetative state[END_REF][START_REF] Kobylarz | Functional imaging of severely brain-injured patients: progress, challenges, and limitations[END_REF][START_REF] Laureys | Auditory processing in the vegetative state[END_REF][START_REF] Laureys | Cortical processing of noxious somatosensory stimuli in the persistent vegetative state[END_REF].

Magnetic Resonance Imaging (MRI) provides information about cerebral morphological lesions as well as information about cerebral functions. Morphological MRI is widely used in patients with chronic disorders of consciousness to evaluate brain injuries. fMRI is based on the measurement of signal modifications induced by the variation of deoxyhemoglobin concentration in the cerebral vasculature (BOLD effect) subsequent to cortical activation [START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF].

It can be performed concomitantly with morphological MRI, within a reasonable time duration. This non-invasive technique does not require any contrast agent injection and it provides localization information. Finally, MR scanners are now available in more and more settings. Although there are technical and ethical limitations in the use of fMRI in the evaluation of patients with chronic disorders of consciousness [START_REF] Kobylarz | Functional imaging of severely brain-injured patients: progress, challenges, and limitations[END_REF][START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF], it could be interesting to evaluate its contribution to the management of patients with acute severe disorders of consciousness.

The aim of this preliminary work was to evaluate the feasability and the potential usefulness of a fMRI approach, when compared to a multimodal approach (EP and PET examinations) in critical care patients.

MATERIAL AND METHODS : PATIENTS :

After ethical committee approval (Comité de Protection des Personnes, Est-IV, Strasbourg), 7 patients (mean age : 49 years (23-73), 3 males, 4 females) presenting with chronic disorders of consciousness in relation with severe brain injuries (prolonged anoxia following out-ofhospital cardiac arrest : 5 patients, severe head traumatism : 1 patient, intracerebral hemorrhage : 1 patient) were included after informed consent obtained from next of kin.

Patients consciousness was evaluated with the Wessex Head Injury Matrix (WHIM) [START_REF] Majerus | Behavioral evaluation of consciousness in severe brain damage[END_REF][START_REF] Laureys | Cortical processing of noxious somatosensory stimuli in the persistent vegetative state[END_REF] and the presence of susceptibility effects in relation to focal brain haemorrhage was evaluated on T2 GE images in the areas elected for multimodal investigation (parietal, temporal and occipital lobes) (see table 1). Patients were free of centrally acting drugs.

METHODS :

The different tests (fMRI, PET and EP) were performed in each patient within a 3-day period of time in a range of 120 (12-554) days after initial brain injury.

MRI :

MRI examination was performed on a 1.5 Tesla imager (Avanto, Siemens, Erlangen, Germany) with an 8-channel head coil. Stimuli were presented using a system dedicated to fMRI (IFIS-SA, Invivo, Orlando, USA). Patients were equipped with headphones for auditory stimuli delivery while protecting from the scanner's noise. A LCD screen was positioned over the head coil in front of the eyes for visual stimuli presentations. In order to limit patients' movements during the examination, their head was fixed using small cushions. Anatomical images were first acquired using an axial T1 weighted 3D GE imaging (TR =11ms, TE =4.94ms, Flip angle =15°, 1 mm 3 isotropic spatial resolution).

Three fMRI runs were acquired using a classical T2*-weighted EPI GE sequence (TR : 4000 ms, TE : 50 ms, Flip angle : 90°, slice thickness : 4 mm, FOV : 256 mm, matrix : 64x64, 120 volume acquisitions/run) with a whole cortex coverage. For each functional run a different perceptual system (somesthesic, auditive and, when possible, visual) was passively stimulated When an activated cluster was identified in one of these expected areas, mean time course and percentage of signal change were further derived in regions of interest (ROIs) located in these areas (defined using boxes of 12mm width) and visually inspected to either confirm or infirm the activation.

Results of fMRI were considered as normal (1) in case of cortical activation in the expected cortical area with time course of signal modification in accordance with variation of BOLD signal, normal with limitation (1*) in case of cortical activation in the expected cortical area but without time course of signal modification in accordance with variation of BOLD signal, absent (0), and not interpretable in case of movement artefacts (NI).

PET :

PET measurement was performed on a PET/CT system (Discovery ST, GEMS, Milwaukee, USA) using the most sensitive 3D mode with in-plane and axial resolutions of 3.27 mm. The CT scanner portion of the Discovery ST consists of a 8-detectors helical CT. Imaging parameters were as follows: 140 kVp, 0 mA, 0.8 s per CT rotation, pitch of 1.675 : 1 mm/s table speed.

Patients were injected intravenously 150 MBq 18 F-FDG at a resting state in a dimly lighted room with minimal background noise. Data acquisition was started 30 min after injection.

Scan duration was 15 min. PET images were reconstructed using CT for attenuation correction with the OSEM algorithm.

PET images were analysed by two readers (CNO, IJN) after coregistration of PET color maps with the 3DT1 (3D Siemens software), in the different functional areas (temporal, occipital and parietal) and cortical fixation was classified as normal (1), hypofixation (0), in case of diminished fixation, or not interpretable (NI) in case of movement artefacts.

EVOKED POTENTIALS :

Patients were tested for 3-modality exogenous evoked potentials (EPs), including visual (VEP), somatosensory EPs (SEP) and event related potentials (ERPs) during a classic auditory oddball paradigm. EPs were recorded with a Nihon Kohden Neuropack8 apparatus (Japan), using standard techniques. All recordings were performed twice to check waveform reproducibility.

SEPs were obtained after average of the response to 200 median nerve electrical stimulations at the wrist with recording sites at Erb's point (to ensure peripheral passage) and above controlateral parietal region. Waves analysis was expressed in terms of presence, deterioration or absence of primary cortical response N20.

VEPs were generated by averaging the responses to 100 reversals of red-light emitting diodes (goggles) presented separately to each eye at a 0.9 Hz frequency, on the recording site located at the median occipital electrode Oz (reference Cz). Wave analysis was expressed in terms of presence, deterioration or absence of occipital primary response.

ERPs were recorded using an auditory oddball paradigm in passive conditions. Deviant (2000 Hz tones; probability of occurrence=0.20) and standard stimuli (1000 Hz tones; probability of occurrence=0.80) were presented binaurally at random with a level of 90 dB. About 50 deviant stimuli were required for each sequence. Responses were recorded with middle frontal Fz-, central Cz-and parietal Pz-electrodes referred to linked earlobes on a 1 sec-base duration of time. Ocular movements or blinking during the session required manually holding of the eyes to decrease artefacts. Wave analysis focused on the presence, the deterioration or the absence of auditory cortical response N100 which corresponds to primary/secondary auditory cortex, as well as of P300 which corresponds to the involvement of more cognitive functions (positive peak occurring in a latency range of 250-600msec in response to deviant stimuli).

Evoked potentials were analysed as normal (1), altered (1*) in case of deterioration of waves, absent (0) or not interpretable in case of movement artefacts (NI).

EVALUATION OF CONCORDANCE BETWEEN fMRI RESULTS, PET AND EP RESULTS:

The results of fMRI activations in somesthesic, visual and auditive cortical areas were compared to the results of EP and to the presence/absence of a normal metabolic activity on PET examination in the same anatomical areas. They were considered as concordant in case of similar results in the same cortical area between fMRI and another examination, non concordant in case of discrepancies between fMRI and another examination or not possible if a result was missing. Because EP have no localizing value (the integrity of one cortical area gives a normal result), we considered as concordant unilateral or bilateral cortical activation observed on fMRI. It was also difficult to distinguish primary and secondary auditive cortex on PET metabolic maps because of poor spatial resolution. This is why we considered as normal the presence of any cortical activation in the auditory cortex (primary or secondary) and compared it to the global metabolic activity in the temporal area.

RESULTS:

Results are presented in Table 2 for our seven patients. Concordances / discordances are recorded as colored boxes in the legend of table 2.

FUNCTIONAL MRI :

fMRI was not interpretable in 4 data-points (patient 6) because of important motion-related artefacts. fMRI was not available in 3 data-points (patients 1,5 and 7), because the patient closed spontaneously their eyes during the examination and this did not allow for visual stimulation. cortical area but without time course of signal modification in accordance with a variation of the BOLD signal.

In seven regions of multimodal analysis (patients 2, 6 and 7) there were susceptibility artefacts observed in relation to cerebral hemorrhage. Among these regions, 4 did not present fMRI activations.

EVOKED POTENTIALS :

EP were not available in one data-point (patient 5 : SEP) and not interpretable in 5 datapoints.

PET :

PET was not interpretable in 3 data-points (patient 4) because anesthesia became mandatory during the examination.

CONCORDANCE BETWEEN FMRI AND OTHER EXAMINATIONS :

Concordance between fMRI and PET: PET results were not interpretable in patient 4, and fMRI results were not either available or interpretable in 7 data-points. This did not allow for comparison of fMRI and PET results in these 11 data-points.

Results were concordant between fMRI and PET in 10 data-points and discordant in 1 datapoint.

Concordance between fMRI and evoked potentials:

EP were not either available in one data-point (patient 5 : SEP) or interpretable in 5 datapoints, and fMRI results were not available or not interpretable in 7 data-points. This did not allow comparison of fMRI and EP results in 13 cases.

Results were concordant between fMRI and EP in 11 data-points and discordant in 4.

Figure 1 displays a concordant pattern of responses (patient #1, upper row) and a discordant pattern of responses (patient #3, lower row).

DISCUSSION

Our results show in many patients a concordance between fMRI and brain functions suggested by EP and resting metabolic activity demonstrated with PET. Moreover, this study demonstrates the feasibility of fMRI in critical care patients, as the number of examinations that are not interpretable is not higher with fMRI (6 data-points), than with either PET (3 datapoints) or EP (6 data-points).

The usefulness of EP and PET scan in the evaluation of brain functions and cerebral metabolism in chronic disorders of consciousness and their potential prognostic value have already been proposed in the literature [START_REF] Laureys | Brain function in coma, vegetative state, and related disorders[END_REF][START_REF] Schiff | Multimodal neuroimaging approaches to disorders of consciousness[END_REF][START_REF] Laureys | Quantifying consciousness[END_REF][START_REF] Owen | Detecting awareness in the vegetative state[END_REF][START_REF] Schiff | Disorders of consciousness[END_REF]. But, although number of recently published studies have evaluated fMRI in this indication, using nociceptive, auditive, visual or more elaborated stimulations [START_REF] Bernat | Can functional MRI detect awareness when a neurological examination does not?[END_REF][START_REF] Giacino | Functional neuroimaging applications for assessment and rehabilitation planning in patients with disorders of consciousness[END_REF][START_REF] Laureys | Brain function in coma, vegetative state, and related disorders[END_REF][START_REF] Laureys | Quantifying consciousness[END_REF][START_REF] Owen | Detecting awareness in the vegetative state[END_REF][START_REF] Boly | When thoughts become action: an fMRI paradigm to study volitional brain activity in non-communicative brain injured patients[END_REF][START_REF] Coleman | Do vegetative patients retain aspects of language comprehension? Evidence from fMRI[END_REF][START_REF] Di | Cerebral response to patient's own name in the vegetative and minimally conscious states[END_REF][START_REF] Metting | Structural and functional neuroimaging in mild-to-moderate head injury[END_REF][START_REF] Owen | Using functional magnetic resonance imaging to detect covert awareness in the vegetative state[END_REF][START_REF] Owen | Residual auditory function in persistent vegetative state: a combined PET and fMRI study[END_REF][START_REF] Perrin | Brain response to one's own name in vegetative state, minimally conscious state, and locked-in syndrome[END_REF][START_REF] Schiff | fMRI reveals large-scale network activation in minimally conscious patients[END_REF], its usefullness as a clinically "easy to perform" examination combining multiple stimulations in association to these gold standards has never been done before in the setting of severe brain injuries. Moreover, our study was prospective, in consecutive non-selected patients, conditions that sometimes lead to the conclusion that a newly available technique should be restricted to specific patients. We report that even in variable and sometimes short period of time between cerebral injury and testing, fMRI evaluation can afford to bring information on remaining brain functioning that is not accessible to clinical evaluation. Although we did not perform all the fMRI in the same range of time after initial damage, this must not be considered as a pitfall of the study, but should rather indicate that this examination has not necessarily to be done as late as in previous studies [START_REF] Owen | Detecting awareness in the vegetative state[END_REF][START_REF] Coleman | Do vegetative patients retain aspects of language comprehension? Evidence from fMRI[END_REF][START_REF] Owen | Residual auditory function in persistent vegetative state: a combined PET and fMRI study[END_REF][START_REF] Schiff | fMRI reveals large-scale network activation in minimally conscious patients[END_REF].

Discordances between fMRI and EP or PET, particularly the absence of fMRI activation in case of normal metabolic cortical activity or normal EP could be related to a different hemodynamic response in comatose patients, when compared to normal subjects, or to the presence of cortical susceptibility effects induced by the presence of hemorrhage.

Our fMRI results can be separated into 3 types : normal cortical activation as compared to those of healthy volunteers, absence of cortical activation, and a group of intermediate results with difficult interpretation because of cortical activation in the expected cortical area but without time course of signal modification in accordance with variation of BOLD signal. In such cases, results need comparison to clinical data, but also to EP and cortical metabolism as evaluated with PET. The discordances between the observed cortical activations and the absence of signal time course modification in relation to the variation of BOLD signal point out one of the main limit of fMRI technique in chronic disorders of consciousness. Indeed the BOLD response is supposed to be based on a normal neurovascular coupling between neuronal activation and local hemodynamics [START_REF] Seiyama | Circulatory basis of fMRI signals: relationship between changes in the hemodynamic parameters and BOLD signal intensity[END_REF]. Vegetative state patients' cortical hemodynamic response could be different, when compared to normal subjects, and could induce aberrant cortical activations [START_REF] Kobylarz | Functional imaging of severely brain-injured patients: progress, challenges, and limitations[END_REF][START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF]. The other limitation of functional MRI is the presence of cortical susceptibility effects, in relation to the presence of cerebral hemorrhage, which could dramatically alter BOLD images. In seven regions of multimodal analysis (patients 2, 6 and 7) there were susceptibility artefacts observed in relation to cerebral hemorrhage. Among these regions, 4 did not present fMRI activations.

The realisation and interpretation of fMRI requires a perfect immobility of the patients (patient 6), but also the absence of centrally acting drugs. This is also the case for PET and EP which are not always easy to perform or to interpret in patients with chronic disorders of consciousness, particularly when anesthesia or analgesia are required (patient 4) for PET or when motion artefacts or electrical interferences occur during the performance of evoked potentials at bedside. Our data suggest that the evaluation of brain functions after severe brain injury would benefit of fMRI because this technique can sometimes compensate for the pitfalls of evoked potentials or PET as indicated in patient 1, in whom auditive and somesthesic EP were not interpretable. Moreover the discordant comparisons between fMRI and PET or fMRI and evoked potentials, are also associated to discordant results between PET and EP (patient 3: vision and auditive). This illustrates the difficulties of evaluation of brain functions in some brain injured patients only with EP and PET: a multimodal approach associating PET, evoked potential and fMRI could therefore be of interest for physicians, patients and even families who desperately expect even tiny information on prognosis . However, MRI is easy to realise, non invasive and supplies in the same time information about morphological brain lesions induced by a trauma or anoxia; moreover, it brings information about brain functions like cortical integration of vision, of audition and of sensitivity. As compared to EP, functional MRI gives a direct visualisation of cortical activation, with good anatomical localisation in both cerebral hemispheres. This information has very important clinical implications since preserved brain functions have a prognostic value [START_REF] Hildebrandt | Brain perfusion and VEP reactivity in occipital and parietal areas are associated to recovery from hypoxic vegetative state[END_REF][START_REF] Naccache | Auditory mismatch negativity is a good predictor of awakening in comatose patients: a fast and reliable procedure[END_REF][START_REF] Giacino | Functional neuroimaging applications for assessment and rehabilitation planning in patients with disorders of consciousness[END_REF]. Indeed knowledge of preserved brain functions could lead to different rehabilitation planning with a better adaptation of stimulation treatment [START_REF] Laureys | How should functional imaging of patients with disorders of consciousness contribute to their clinical rehabilitation needs?[END_REF][START_REF] Lew | Prognostic value of evoked and event-related potentials in moderate to severe brain injury[END_REF]. For example, patient 1 had no visual response but some auditory integration response; therefore clinical rehabilitation was axed on auditory stimulation. These data are important for the medical staff in charge of the patients.

In conclusion, fMRI increases the accuracy of neurological evaluation in our series of critically ill patients with acute brain injury, provided that the BOLD signal can be properly interpreted. It turns out, to be an additional tool of interest to complete evaluation in those patients where EP and PET examinations are discordant. In addition it may help physicians to evaluate with more certainty severe residual neurological impairment and progress in ethical choices if a real prognostic value can be observed among large cohorts of patients.

  following a block paradigm design. Block paradigms consisted in 6 alternations of 10 volumes either with or without stimulation. Somesthesic stimulation consisted in passive tactile mechanical brushing of alternatively right and left legs. Auditive stimulation was performed by the binaural delivery of an 2000 Hz pure tone pulsed at 1Hz through the headphones. When patients had spontaneously their eyes open (4 subjects), visual stimulation was performed through presentation of a 8 Hz reversing checkerboard on the LCD screen. Stimuli presentation was semi-randomly ordered among patients. Functional images were individually processed using SPM2 software (Wellcome Department of Cognitive Neurology, Queen's College, London; http://www.fil.ion.ucl.ac.uk/spm/). The pre processing included a rigid body motion correction, a co registration on the anatomical image using mutual information method, and smoothing by a 8mm isotropic Gaussian kernel. Examinations with more than 2 mm translation or 6 ° rotation movements were excluded from further analysis. Individual statistical analyses were performed by solving the linear generalized model in which regressors of interest were modelled by a boxcar function representing the time-course of the stimuli delivery convolved with the canonical hemodynamic response function provided by SPM. Motion parameters were included in the model as regressors of no interest. A high pass filter allowed for the removal of slow drifts in the baseline. Student t-test were computed voxelwise for each stimulus of interest and thresholded at p<0.005 providing a stimulus related activation map. Obtained cortical activation maps were then fused with each patient volumetric T1 anatomical image. Cortical activation areas of somesthesic, auditive and visual stimulations are well established (27-29). In healthy subjects, somesthesic stimulation of the leg induces activation of primary sensitive cortex situated in the medial part of contralateral postcentral gyrus (Broadman area (BA) : 1, 2 and 3 ). Bilateral auditive stimulation induces activation of bilateral primary auditive cortex situated in Heshl gyrus (BA 41, 42) and secondary auditive cortex situated around of it (BA 22). Bilateral visual stimulation induces activation of bilateral primary visual cortex situated along the calcarine sulcus and in the inferior occipital gyrus (BA 17, 18 and 19).

Figure 1 :

 1 Figure 1 : Patient 1 : Figure 1a (Auditive Functional MRI), 1b ( 18 F-FDG PET scan), 1c

Table 1 :

 1 Patients clinical and radiological characteristics. CA : Cardiac Arrest, HT : Head Traumatism, IH : Intracerebral Hemorrhage, MCS : Minimally Conscious State, VS : Vegetative State, Time: Time from acute injury to test performance (in days), cerebral hemorrhage : presence of susceptibility effects in relation to intracerebral hemorrhage on T2 GE sequence.

		Patient 1	Patient 2		Patient 3	Patient 4	Patient 5	Patient 6	Patient 7
	Age (years)	35	63		67	33	73	52	23
	Etiology	CA	CA		CA	CA	CA	HT	IH
	Clinical								VS
	evaluation	of
	consciousness	VS	MCS		MCS	MCS	VS	VS
	WHIM	0	27		21	22	0	7	2
	Time (days)	12	32		143	14	20	554	68
			Yes : right and					Yes : right
			left occipital				Yes : left	temporal and
	Cerebral		lobe,	left				temporal and	occipital lobe
	hemorrhage	no	parietal lobe	no	no	no	parietal lobe

Table 2 :

 2 Synthesis of results :NA : not available. PET : 0=hypofixation, 1=normal, NI= not interpretable. Evoked potentials : SEP : somatosensory evoked potentials, VEP : visual evoked potentials, 1=normal, 1*=altered, 0=absent, NI : not interpretable. fMRI : functional MRI, fMRI I : primary auditive area, fMRI II : secondary auditive area, 1=normal, 0=Absent, 1*= normal with limitation.

Somatosensor y

Conflict of interest : None


Keywords: Lamb waves, Reflection, Viscosity, Viscous fluids

This paper is a study of a new sensor for fluid characterization. This sensor is composed of a stainless steel plate in contact with a viscous material. The aim is to characterize the material viscosity by using reflected Lamb waves at the boundary interface. In order to identify the effects on the Lamb reflected modes by the viscous material, a complete study of the propagation wave in the alone plate is first presented. The propagation modes of the loaded plate are then investigated. By monitoring the mechanical impedance, the viscosity of the material in contact is extracted. In order to validate the experimental set-up, the mechanical impedance variation is measured for different water-glycerol mixtures. Results are in good agreement with those obtained by other techniques in the literature.

I. INTRODUCTION

wide range of commercial products have a great interest for their mechanical, electrical, optical and biocompatibility properties because of their colloidal structures. These include food products (such as milk, yogurt, mayonnaise), sprays and pharmaceutical products. These products exist either in sol phase or pass through a colloidal state during their manufacture. The monitoring of these complex fluids evolution is singularly crucial during the first step of their formation. Several techniques can be used to characterize the structure evolution [START_REF] Yan | In-situ determination of particle size distributions in colloids[END_REF]- [START_REF] Bike | Measuring colloidal forces using evanescent wave scattering[END_REF]. To monitor the mechanical properties of these materials, low frequency rheometers are commonly used. They give access to the complex dynamic shear modulus [START_REF] Ponton | Rheology of titanium oxide based gels: determination of gelation time versus temperature[END_REF]- [START_REF] Griesmar | Kinetic study of silicon alkoxides gelation by acoustic and rheology investigations[END_REF] during the transition. More recently acoustic investigation have been carried out for the study of the sol-gel materials evolution [START_REF] Senouci | Acoustic resonance in tetramethoxysilane matrices: A new tool to characterize the gel formation[END_REF]. Due to the frequency range used by these techniques, the first step of the microscopic evolution cannot be detected. Ultrasonic methods based on the measurement of either the propagation speed of acoustic waves or their attenuation by ultrasonic spectroscopy [START_REF] Matsukawa | An ultrasonic evaluation of the sol-gel synthesis of silica glass[END_REF] can also be performed. Theoretical works have been Manuscript received February 11, 2009. N. Wilkie-Chancellier, L. Martinez, S. Serfaty and P. Griesmar are with the University of Cergy-Pontoise, ENSEA, Equipe Circuit Instrumentation et Modélisation en Electronique (ECIME-EA04), F-95000 Cergy-Pontoise, France. (e-mail: Nicolas.Wilkie-Chancellier@u-cergy.fr). performed to study the interaction of waves with a viscous fluid [START_REF] Zhu | The propagation of Lamb waves in a plate bordered with a viscous liquid[END_REF]- [START_REF] Bianco | Influence of viscosity on the reflection coefficient at a fluid-solid interface[END_REF]. On the other hand, measurements of the viscosity [START_REF] Shah | Effect of viscosity on ultrasound wave reflection from a solid/liquid interface[END_REF] or the viscosity-density [START_REF] Lévêque | An acoustic sensor for simultaneous density and viscosity measurements in liquids[END_REF]- [START_REF] Greenwood | Measurement of the viscosity-density product using multiple reflections of ultrasonic shear horizontal waves[END_REF] product have been carried out using ultrasonic shear wave [START_REF] Shah | Mesuring Newtonian viscosity from the phase of reflected ultrasonic shear wave[END_REF]- [START_REF] Greenwood | Measurement of viscosity and shear wave velocity of a liquid or slurry for on-line process control[END_REF]. Ultrasonic guided waves have been used for viscosity measurements in melt using in a alumina buffer rod [START_REF] Prasad | Viscosity measurements of melts at high temperatures using ultrasonic guided waves[END_REF].

Works have been led on the study of viscoelastic anisotropic materials [START_REF] Hosten | FE modeling of Lamb mode diffraction by defects in anisotropic viscoelastic plates[END_REF]- [START_REF] Hosten | Comments on the ultrasonic estimation of the viscoelastic properties of anisotropic materials[END_REF]. Recently, other studies based on the resonance of an AT cut quartz sensor have been developed. Indeed, the characterization of complex viscoelastic fluids is possible using shear wave propagation in the complex fluid [START_REF] Cernosek | Comparison of lumped-element transmissions line models for thickness-shear-mode quartz resonator sensors[END_REF]- [START_REF] Bandey | Modelling the response of thickness-shear-mode resonators under various loading conditions[END_REF]. Therefore a complete monitoring has been carried out from liquid to gel state. This experimental technique gives information from the microscopic to the macroscopic scale. For example, restructuration of yogurt (i.e. casein network) has been pointed out just before gelation [START_REF] Farrell | Particles sizes of casein submicelles and purified k-casein -Comparisons of dynamic light scattering and electron microscopy with predictive three-dimensional molecular models[END_REF]- [START_REF] Ould-Ehssein | Ultrasonic monitoring of yoghurt formation by using AT-cut quartz: lighting of casein micelles interactions process during the acidification[END_REF]. Moreover, other works have clearly shown the interest of a new viscoelastic time (precursor to the gelation time) characteristic of sol-gel material evolution [START_REF] Ould-Ehssein | Kinetic study of silica gels by a new rheological ultrasonic investigation[END_REF].

However, for all these techniques based on bulk waves propagation, the sensor has to be placed into the studied material. In order to deport the sensor, it seems to be interesting to use the properties of the surface wave propagation. Recently, applications have been studied by the mean of a microfluidic surface acoustic wave sensor platform based on Love wave propagation [START_REF] Raimbault | A microfluidic surface acoustic wave sensor platform : Application to high viscosity measurements[END_REF]. But up to now experimental techniques based on the Lamb wave propagation are not widespread in this subject. However, the reflection of a Lamb wave at the free edge of a straight cut plate has been theoretically

[23]- [START_REF] Rose | Ultrasonic waves in solid media[END_REF], numerically [START_REF] Leclezio | Numerical predictions and experiments on the free-plate edge mode[END_REF]- [START_REF] Wilkie-Chancellier | Numerical description of the edge mode at the beveled extremity of a plate[END_REF] and experimentally [START_REF] Wilkie-Chancellier | Experimental study and signal analysis in the Lamb wave conversion at a bevelled edge[END_REF] investigated.

The aim of this paper is to develop a sensor and an experimental technique based on the Lamb wave reflection phenomena in order to characterize viscoelastic materials. The great interest of this sensor is to assume a remote viscoelastic measurement. This presented work is a preliminary study on Newtonian liquid to validate the experimental device and the associated model.

In a first time, the reflection phenomena are experimentally investigated if a Lamb mode is incident on a free end. The propagation of the Lamb waves is followed using of a laser vibrometer. At the end of the plate, several refleted modes are qualitatively and quantitatively determinated from a 2D Fast Fourier Transforms signal treatment.

In second step, an experimental study is carried out at the

Lamb Wave Sensor For Viscous Fluids Characterization

Nicolas Wilkie-Chancellier, Loïc Martinez, Stéphane Serfaty, Member IEEE, and Pascal Griesmar A end of the plate in contact with several water-glycerol mixtures. Referencing to the study at the free edge, the acoustic impedance of the viscous fluid is computed using the measured reflection coefficient at the plate-mixture interface.

Experimental results are compared with the theoretical Newtonian liquid model. This model is suitable for the characterization of viscoelastic materials.

II. LAMB WAVE SENSOR INVESTIGATION

A. Sensor Description Lamb waves are generated using a wedge transducer on a stainless steel plate (density ρ =7800 kg/m 3 , longitudinal and transversal velocities respectively c L = 5850 m/s and c T = 3150 m/s). This plate is 150 mm long, 30 mm wide and 2 mm thick (thickness E). Its extremity is a normal free edge (Fig. 1). The transducer is made of a 1.5 mm thick piezoelectric plate (25x35 mm) and a plexiglass wedge to have a θ angle with the plate. The wedge is placed at d =10 cm from the plate end. This wedge transducer is chosen to generate a linear Lamb wave front (25 mm wide) at the resonant frequency F = 1.25 MHz. The frequency-thickness product FE is then equal to 2.5 MHz mm. At this frequency-thickness product, two different antisymmetric modes (A 0 and A 1 ) and one symmetric mode (S 0 ) are possible. In addition with the normal free edge, there is no change of symmetry modes at the reflection.

Taking into account the plexiglass properties, the wedge angle is chosen to 24 degrees to generate the A 1 Lamb mode because this mode is more sensitive to the conversion rate.

For our application, this sensor is used in reflectometry. Then, the reflected A 1 Lamb mode can only be received.

B. Lamb Wave Propagation in the Plate for Free Boundary

To characterize a fluid in contact at the end of the plate, the understanding of the unloaded sensor behavior is required. In fact, the knowledge of the reflection phenomena gives access to the interactions of surface waves with the fluid to characterize. To achieve this stage, the surface wave propagation is studied.

In order to study the propagation of Lamb waves in the plate, the experimental set-up shown in Fig. 2 is implemented. Three units can be distinguished. The first one generates the A 1 Lamb mode using our transducer design. In order to generate this wave, a burst signal of twenty periods is applied to the transducer. A Polytec He-Ne laser vibrometer (OFV 505) coupling with a micro-positioning system scans by interferometry the normal surface plate displacements U 2 along the propagation. For each position of the vibrometer, its time response is acquired using a Lecroy digital scope (WS 424) and stored in a computer. Along the last 5 cm of the plate the U 2 displacements are collected every 0.2 mm step. The Fig. 3 shows the spatiotemporal signal s(x 1 ,t) of U 2 . The signal observed on the figure results from the superposition of the incident wave and the reflected wave as:

( ) ( ) ( ) x k t j r x k t j i e U e U t x U + - + = ω ω , 2 (1.)
where U i and U r . are respectively the incident and reflected wave magnitudes and k is the complex wave number.

Note that the reflection coefficient R can be extract from (1) and can be written by:

x k j i r e U U R 2 = ( 2.
)

In order to determinate the reflected waves, two successive Fourier transforms (temporal and spatial) of the total signal s(x 1 ,t) are computed to obtain the Ψ(k,FE) signal [START_REF] Alleyne | A two-dimensionnal Fourier transform method for the measurement of propagating multimode signals[END_REF]. The incident and reflected Lamb waves can be observed in the dual space, respectively for k>0 and k<0 in Fig. 4. The theoretical dispersion curves of the Lamb waves have been reported in dashed lines on this figure. The different reflected modes (A 1 and the A 0 Lamb modes) can be then distinguished. The quantitative study is performed by measuring each magnitude mode and by computing energy balance for the given incident A 1 Lamb mode [START_REF] Morvan | Lamb wave reflection at the free edge of a plate[END_REF]. At FE = 2.5 MHz mm, 77% of the total energy is converted into A 0 mode and 23% into A 1 .

C. A 1 Lamb Mode Evolution Versus Fluid Loading

If the edge of the stainless steel plate is in contact with a fluid, the energy balance is modified. According to the boundary conditions at the plate-fluid interface, the reflection coefficient magnitude depends of the fluid characteristics [START_REF] Aristégui | Guided waves in fluid-filled pipes surrounded by different fluids[END_REF]- [START_REF] Peter | Viscosity-induced attenuation of longitudinal guided waves in fluid-loaded rods[END_REF].

In order to confirm that the A 1 mode is significantly modified by the boundary conditions, the Lamb wave evolution is monitored for water and for glycerol. These fluids are chosen to monitor the A 1 mode evolution in a wide range of absolute viscosity (respectively η = 1.005 mPa.s and η = 1.780 Pa.s for water and glycerol).

Table I shows the modifications of the energy balance of the A 1 reflected Lamb mode. These preliminary results validate that our device can be used to remotely characterize fluids.

III. VISCOUS FLUID CHARACTERIZATION

A. Experimental Set-up for Fluid Characterization

Now, the Lamb wave sensor is used to characterize viscous fluids in contact with the stainless steel plate. The wedge transducer is used in reflectometry. The plate can be then considered as a transmission line. The experimental set-up is shown in Fig. 5. The reflection coefficient induced by the water-glycerol mixture can be extrated from the reflected mode magnitude U r . Indeed, at the plate-fluid interface, the reflection coefficient expression is written as :

d k j i r e U U R 2 = (3.)
where i U is the magnitude of the incident A 1 Lamb mode and d is the distance from the transducer to the plate end. This magnitude U i only depends on the magnitude of the electrical pulse and the mechanical conditions of the contact between the wedge transducer and the steel plate.

Then, the U i can be experimentally deduced from the free edge configuration (when keeping constant the incident amplitude in both experimental studies). The Lamb wave reflection coefficient is written as :

air r r R U U R air = (4.)
where R air is the reflection coefficient for free edge configuration which is expressed as : (5.)

B. Extraction Of The Viscous Fluid Characteristics

In order to characterize the mechanical properties of the fluid in contact with the steel plate, the acoustic impedance Z m of the fluid can be extract from (5) introducing the acoustic impedance Z steel of the stainless steel:

R R Z Z steel m + - = 1 1 (6.
)

Introducing ( 5) in ( 6), the acoustic impedance of the fluid can be expressed as a function of the experimental parameters:

air r r air r r steel m R U U R U U Z Z air air + - = (7.
)

Taking into account that the boundary conditions at the plate-fluid interface generate a shear stress the acoustic impedance Z m of the fluid depends on the complex shear modulus [START_REF] Christensen | Theory of viscoelasticity : an introduction[END_REF]. For Newtonian liquid the acoustic impedance Z m can be then linked to dynamic viscosity η of the fluid and the density ρ of the fluid:

( ) 5 . 0 ω η ρ j Z m = (8.)
where ω is the resonant angular frequency of the wedge transducer. Then, the acoustic impedance of the viscous fluid is proportional of the square root of the ρ η product.

Sensors-3112-2009.R1

IV. RESULTS AND COMPARISON

Different water-glycerol mixtures are used to validate our fluid characterization device. The temperature fluid is maintained at 25°C by a thermostated cell. At this temperature the well known density ρ (in kg/L) and dynamic viscosity η (in mPa.s) are given in Table II.

In this table, the quantity of glycerol in water is given by the weight ratio x m,gly defined as follow: where m gly and m water are respectively the glycerol and water weight. Note that water without glycerol is written x m,gly =0% and glycerol without water x m,gly =100%.

Water-glycerol mixtures used in this study are in weight ratio range from 0% to 100% with 10% step.

To experimentally extract the viscosity of mixtures, the acoustic impedance of fluid Z m is first measured from the evolution of the reflection coefficient. The figure 6 shows the evolution of Z m versus the density viscosity product (ρη) 0.5 . This figure shows that the measurements of absolute viscosity are in good agreement until with the tabulated values from (ρη) 0.5 =2.76 kg/m 2 /s 0.5 (i.e. x m,gly =50%).

From Z m measurements, the experimental density-viscosity product is therefore computed. A comparison with the tabulated (ρη) 0.5 values are reported in Table II. Consequently our Lamb wave sensor can accurately measure the absolute viscosity of weakly viscous fluid (η ≤ 7 mPa.s). This domain of validity is comparable to the usual techniques in the literature, such as using multiple reflections of ultrasonic shear horizontal waves [START_REF] Greenwood | Measurement of the viscosity-density product using multiple reflections of ultrasonic shear horizontal waves[END_REF][START_REF] Shah | Mesuring Newtonian viscosity from the phase of reflected ultrasonic shear wave[END_REF][START_REF] Greenwood | Measurement of viscosity and shear wave velocity of a liquid or slurry for on-line process control[END_REF] or using the phase shift of the reflected ultrasonic shear wave [START_REF] Shah | Mesuring Newtonian viscosity from the phase of reflected ultrasonic shear wave[END_REF]. For more viscous fluids (ie 60% ≤ x m,gly ≤ 100%), the error increases, pointing out the limits of the viscosity measurement with Lamb waves. Because of the dominant nature of the Lamb wave propagation, the determination of the high bulk viscosity becomes difficult due to the high dissipation of the shear wave energy in the fluid. In our model, the steel plate is considered as a lossless transmission line. In fact, at the ε depth contact layer between the stainless steel plate and the mixture (Fig. 1), the effective propagating medium depends on the local interactions (such as wettability , fluid-plate meniscus) and on the mechanical properties of the both materials in contact. Equations ( 3) and (4) should include these interactions which strongly increase for high η values.

Using Table II, such a steep transition of the dynamic viscosity occures for x m,gly ≥ 60%. A possible way to improve the accuracy is to take into account these propagation effects.

V. CONCLUSION

The experimentation performed with the A 1 Lamb mode shows that guided waves are able to detect the change of the viscous fluid parameters from their reflection coefficient. The acoustic impedance of the viscous material is expressed as a function of the reflection coefficient, and more particularly as a function of the experimental parameters. The experimental study with water-glycerol mixtures allows us to determine the absolute viscosity for weak viscous fluids. For Newtonian liquids, the extracted dynamic viscosity is in good agreement with the tabulated data for small values less than 7 mPa.s. The high dissipation of the shear wave energy in the strongly viscous fluids could explain the divergence of measurement.

However, taking into account the complex shear modulus the Lamb wave sensor presented in this paper should be used to characterize the viscoelasticity of fluids. This technique can be then a new way to monitor the first steps of weak gels formation.

Fig. 1 .

 1 Fig. 1. Description of the Lamb wave sensor.

Fig. 2 .

 2 Fig. 2. Experimental set-up for Lamb mode determination (measurements using the vibrometer).

Fig. 3 .

 3 Fig. 3. Time-space representation of the normal displacement U2 on the plate.

Fig. 4 .

 4 Fig. 4. Dual space representation of the normal displacement on the plate.

Fig. 5 .

 5 Fig. 5. Experimental set-up for fluid characterization.

Fig. 6 .

 6 Fig.6. Comparison of the experimental acoustic impedance for waterglycerol mixtures versus tabulated values (the squares are the experimental data points, the continuous curve is the best fit using a logarithm function and the dashed line is the tabulated values).

TABLE II TABULATED

 II VERSUS EXPERIMENTAL (ρη) 0.5 VALUES FOR WATER-GLYCEROL MIXTURES AT 25°C

	Glycerol in water xm,gly	Tabulated ρ (kg / L)	Tabulated η (mPa.s)	Tabulated (ρη) 0.5 ( kg / m 2 / s 0.5 )	Experimental (ρη) 0.5 ( kg / m 2 / s 0.5 )	Error (%)
	0 %	1.000	1.005	1.002	0.954	4.7
	10 %	1.023	1.341	1.171	1.214	3.6
	20 %	1.048	1.844	1.390	1.376	1.0
	30 %	1.073	2.683	1.697	1.745	2.8
	40 %	1.100	4.093	2.122	2.199	3.6
	50 %	1.127	6.762	2.761	2.676	3.1
	60 %	1.155	12.474	3.796	3.369	11.2
	70 %	1.182	26.595	5.607	4.062	27.6
	80 %	1.210	72.721	9.380	4.415	52.9
	90 %	1.238	271.122	18.321	6.106	66.7
	100 %	1.263	1780.830	47.426	7.506	84.2

Nicolas Wilkie-Chancellier was born on July 26, 1976, in Le Havre, France. He studied physics, acoustics and their applications at the University of Le Havre, France. He obtained the Ph.D. degree in 2003 on the topic "Reflection and conversion of a Lamb wave at a bevelled end of a plate" in the University of Le Havre. He continued his research on the health monitoring applications in aeronautics, cultural heritage and civil engineering using Nonlinear Elastic Wave Spectroscopy as a post-doctoral researcher in Katholieke Universiteit Leuven, campus of Kortrijk, Belgium. Assistant Professor in the University of Cergy-Pontoise since 2004, his research activities are focused on NDT using surface and bulk waves, characterization of complex materials, and acoustic instrumental systems design. He is member of the GAPSUS group of the French Acoustic Society (SFA).

Loïc Martinez was born in 1967 in


Adjih -E Baccelli -P

-P Jacquet

-M Minet

Philipp -G

Wittenburg

Février

C Adjih

E Baccelli

P Jacquet

P Minet

M Philipp

G Wittenburg

INFORMATIQUE ET EN AUTOMATIQUE

Keywords: Sensor, Ad hoc, Wireless, Network, IETF, Internet, Things, Routing, Deployment

Protocols that are to be employed in the context of the Internet of Things (IoT) have to meet a wide variety of application-specific requirements [3] [4] [1] [2]. In this report, we reflect on recent experiences, gained from several real-world deployments in which we have participated, which use low power, embedded networking devices. We discuss the lessons learned from these deployments, with an emphasis on questions affecting the IP layer and, in particular, on the routing protocols for these networks. We point out open issues and possible directions of future work for such routing protocols.

Quelques Déploiments de Réseaux de Capteurs

Tests pour les Normes IETF concernant l'Internet des Objets

Résumé : Ce document rapporte certaines observations à propos de déploiments de réseaux de capteurs dans le contexte de l'Internet des Objets, ayant trait au protocole de routage utilisé dans ce contexte, et notamment la norme en la matière, actuellement en cours de développement à l'IETF.

Mots-clés : Capteurs, Internet, Objets, IETF, Routage, Réseau, Déploiment, Ad Hoc

Protocols that are to be employed in the context of the Internet of Things (IoT) have to meet a wide variety of application-specific requirements [START_REF] Dohler | Routing Requirements for Urban Low-Power and Lossy Networks[END_REF]

[4] [1] [2].
In the following, we reflect on recent experiences, gained from several real-world deployments in which we have participated, which use low power, embedded networking devices. We discuss the lessons learned from these deployments, with an emphasis on questions affecting the IP layer and, in particular, on the routing protocols for these networks. We point out open issues and possible directions of future work for such routing protocols.

Deployment Experience

As part of our work on the SensLAB1 project, we participate in a large, general-purpose testbed for wireless sensor networks built by INRIA. SensLAB aims to be an accurate and efficient scientific tool that can be employed for the design, development, tuning, and experimentation of real-world, large-scale sensor network applications. Since 2009, SensLAB has been composed of 1,024 nodes distributed across four sites in France. The SensLAB nodes are based on the Texas Instruments MSP430 micro-controller and employ one of two radio interfaces, one operating in the 868 MHz ISM band and the other at 2.4 GHz, both in compliance with the IEEE 802.15.4 standard. We have conducted experiments on such hardware using various routing protocols including OLSR, and RPL [START_REF] Winter | RPL: IPv6 Routing Protocol for Low power and Lossy Networks[END_REF] (including the P2P extension [START_REF] Goyal | Reactive Discovery of Point-to-Point Routes in Low Power and Lossy Networks[END_REF]) on top of the Contiki operating system [START_REF] Dunkels | Contiki -a Lightweight and Flexible Operating System for Tiny Networked Sensors[END_REF].

For RPL, we found that the software stack, i.e., Contiki, 6LowPAN, IPv6, RPL (upward routes, no DAOs), and the P2P extension, barely fits into the 48 KB of flash memory available on the MSP430 (current state of the art for such devices), leaving close to no space for any application-specific code. Our experiments showed the Trickle mechanism for establishing upward routes works efficiently and converges within reasonable time. However, the way the RPL specification handles asymmetric links (common in wireless networks) required special attention from our part. Depending on the application requirements, relying on NUD may be inappropriate because a path is required at the time data needs to be sent. In these cases, appropriate link metrics must be used in addition to the RPL specification to continuously monitor link availability. In our experiments, some metrics were functional at small scale (such as ETX), but more work is needed to determine which metrics/signalling are appropriate for larger, denser networks, and that fit the memory/energy limitations on each device.

We also took part in is the OCARI (Optimization of Communications in Ad hoc Industrial Networks) 2 project, started in 2006, with partners including international electricity provider EDF, DCNS, and Telit-RF Technologies. This project targets the monitoring of industrial equipment or civil engineering constructions using wireless sensor networks, including performance testing of equipment, radioprotection of site maintenance, and state control of devices. The application requirements in this project include time-bounded delays for specific types of traffic, support for nomadic nodes (moving at pedestrian speed) for data collection, energy efficiency to maximize network lifetime, and ease of deployment. Before ROLL was chartered, this project had us design an energy efficient routing protocol, and a node activity scheduling algorithm based on node coloring to save energy by allowing nodes (including routers) to sleep without incurring any data loss.

In the context of the AVS-Extrem3 project, we have developed a wireless sensor platform for perimeter surveillance. The deployment scenario consists of embedded sensor nodes built into the metallic structure of fence elements. These nodes sample acceleration data, use distributed pattern matching algorithms to detect events, and report security-relevant events to a base station. Network traffic in this scenario consists of rare event notification packets sent from sensor nodes to the base station. Reception of these packets needs to be acknowledged by the base station. Due to the widespread deployment of devices along the perimeter of the area under surveillance, a routing protocol needs to operate in a hop-by-hop fashion, i.e., source routing would result in prohibitively large packet headers.

In the context of the ASIST european project (to start in 2011, pending approval), a large deployment is planned, interconnecting heterogeneous sensor networks over an IP backbone, targeting airport monitoring. Each individual sensor node will have to be accessible from the Internet, and reliable node-to-node communication across network borders will be required. Furthermore, embedded devices will have to support enumeration and service discovery.

Lessons Learned and Position

The deployments described in the previous section are well within the IoT scope. Yet, the requirements they impose on the routing protocol are varied, including some requirements that are not targeted by current RFCs or RFCs-to-be, such as sensor mobility, or traffic patterns requiring paths from the sink to each sensor/actuator in the network. We submit that it is necessary to address these requirements as they are not uncommon.

In order to address these further requirements, we suggest that RPL be complemented with additional extensions. However, we stress that the resourceconstrained nature of the targeted devices is a concern for such extensions-tobe-designed, since our experience shows that the RPL specification is already complex, and that available network stacks including an RPL implementation barely fit in the memory of state-of-the-art sensors. Furthermore, the interoperability of the various existing or future extensions and modes of operation should be addressed explicitly in the RPL specification. Finally, we submit that it is important that the basic RPL specification addresses properly the issue of unidirectional links when NUD is not usable.

http://www.senslab.info/

http://www.ocari.org/ RR n°7551

http://cst.mi.fu-berlin.de/projects/AVS-Extrem/ INRIA


Sourour Ammar
email: sourour.ammar@univ-nantes.fr

Philippe Leray
email: philippe.leray@univ-nantes.fr

François Schnitzler
email: fschnitzler@ulg.ac.be

Louis Wehenkel
email: l.wehenkel@ulg.ac.be

Sub-quadratic Markov tree mixture learning based on randomizations of the Chow-Liu algorithm

The present work analyzes different randomized methods to learn Markov tree mixtures for density estimation in very high-dimensional discrete spaces (very large number n of discrete variables) when the sample size (N ) is very small compared to n. Several subquadratic relaxations of the Chow-Liu algorithm are proposed, weakening its search procedure. We first study naïve randomizations and then gradually increase the deterministic behavior of the algorithms by trying to focus on the most interesting edges, either by retaining the best edges between models, or by inferring promising relationships between variables. We compare these methods to totally random tree generation and randomization based on bootstrap-resampling (bagging), of respectively linear and quadratic complexity. Our results show that randomization becomes increasingly more interesting for smaller N/n ratios, and that methods based on simultaneously discovering and exploiting the problem structure are promising in this context.

Introduction

Directed probabilistic graphical models encode a joint distribution over a set of variables by a product of conditional probability distributions, one for each variable conditionally to its parents in the directed graph. These models may be learned from data and used to perform probabilistic inferences over the encoded distribution [START_REF] Pearl | Fusion, propagation, and structuring in belief networks[END_REF]. However, exact inference and learning with such models are both NPhard, unless the skeleton of the graph is constrained [START_REF] Cooper | The computational complexity of probabilistic inference using bayesian belief networks[END_REF]. Existing learning algorithms are not scalable to high dimensional spaces because of their excessive computational complexity [START_REF] Auvray | On the construction of the inclusion boundary neighbourhood for Markov equivalence classes of Bayesian network structures[END_REF].

Markov Trees are an interesting subclass of directed graphical models, whose skeletons are acyclic and for which each node of the graph has (at most) one parent. With Markov trees, the computational complexity of probabilistic inference and parameter learning are linear in the number of variables [START_REF] Pearl | Fusion, propagation, and structuring in belief networks[END_REF]. Further, Markov tree structures may be learned efficiently by the Chow-Liu algorithm (Section 3.1), quadratic in the number of variables.

While Markov trees impose strong modeling restrictions, mixtures of Markov trees can represent a much wider (actually unlimited) class of probability densities than single Markov trees while retaining their interesting computational properties in terms of inference [START_REF] Meila | Learning with mixtures of trees[END_REF], making these models attractive for scaling graphical models to high-dimensional problems. As a matter of fact, these simple graphical models were used for these reasons, in order to build optimized mixtures of models for probability density estimation [START_REF] Meila | Learning with mixtures of trees[END_REF] by using an expectationminimization algorithm.

In supervised learning, a generic framework which has led to many fruitful innovations is called "Perturb and Combine". Its main idea is to on the one hand perturb in different ways the optimization algorithm used to derive a predictor from a dataset and on the other hand to combine in some appropriate fashion an ensemble of predictors obtained by multiple iterations of the so perturbed search algorithm. This approach may in particular lead to a strong reduction in variance [START_REF] Breiman | Arcing classifiers[END_REF]. It was first explored for probability density estimation in [START_REF] Ammar | High-dimensional probability density estimation with randomized ensembles of tree structured Bayesian networks[END_REF] by comparing various kinds of large ensembles of simple graphical models in the form of Markov Trees.

The above mentioned algorithms for learning mixtures of Markov trees use the Chow-Liu algorithm [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]. However, since this algorithm is quadratic in the number of variables, these methods do not scale well to very high-dimensional problems, with thousands or even millions of variables. Thus, [START_REF] Schnitzler | Towards sub-quadratic learning of probability density models in the form of mixtures of trees[END_REF][START_REF] Ammar | Subquadratic Markov tree mixture models for probability density estimation[END_REF] tried to investigate mixtures of models learned by using various randomized versions of the Chow-Liu algorithm, with the aim of reducing the computational complexity below the quadratic level, and simultaneously improving accuracy in small (i.e. realistic) sample size conditions by variance reduction. The aim of the present paper is to analyse these methods and compare their results in the same framework.

The rest of the paper is organized as follows. We describe tree models and models of mixtures of trees more formally in section 2, and state in section 3 the different tree mixture learning algorithms that we want to analyse. We explain and discuss our empirical evaluation of these algorithms in section 4, before concluding.

Mixtures of Markov trees

Let X = {X 1 , . . . , X n } be a finite set of discrete random variables, and

D = (x 1 , • • • , x N )
be a sample (we will use the term "dataset" to denote it) of joint observations

x i = {x i 1 , • • • , x i n }
independently drawn from some data-generating density P G (X 1 , . . . , X n ).

A mixture distribution P T (X 1 , . . . , X n ) induced by a multiset T = {T 1 , . . . , T m } of m Markov trees is defined as a convex combination of elementary Markov tree densities, i.e.

P T (X) = m i=1 µ i P T i (X),
where µ i ∈ [0, 1], m i=1 µ i = 1, and P T i (X) is the probability density over X encoded by the graphical model composed of the Markov tree structure S i and its parameter set θi :

P T i (X) = P S i , θi (X) = n p=1 P θi (X p |P a S i (X p )),
where P a S i (X p ) is the parent variable of X p in the tree structure S i .

Several versions of Markov tree mixtures were studied in [START_REF] Ammar | Probability density estimation by perturbing and combining tree structured Markov networks[END_REF][START_REF] Ammar | High-dimensional probability density estimation with randomized ensembles of tree structured Bayesian networks[END_REF] as an alternative to classical methods of density estimation in the context of highdimensional spaces and small datasets: mixtures of tree structures generated in a totally randomized fashion with linear complexity in the number of variables and ensembles of optimal trees derived from bootstrap replicas of the dataset by the Chow and Liu algorithm [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]) (i.e. bagging of Markov trees).

Other studies tried to relax the Chow-Liu algorithm to reduce its computational complexity while maintaining its accuracy [START_REF] Schnitzler | Towards sub-quadratic learning of probability density models in the form of mixtures of trees[END_REF][START_REF] Ammar | Subquadratic Markov tree mixture models for probability density estimation[END_REF]. In the present work we analyze these methods in terms of computational complexity, accuracy, and running time, within a same framework.

Panel of learning algorithms

Algorithm 1 describes our general methodology to learn a mixture of m Markov trees from a dataset D.

It may be declined by using different variants of the three subroutines it uses, namely Build-MarkovTreeStructure, LearnPars, and Comp-Weights. In this paper we focus on the effect of varying only the first one, used to infer the structures S i of the mixture terms. Next, we describe the different versions of Build-MarkovTreeStructure that we have considered in our study. We start by describing the original Chow-Liu method.

Algorithm 1 (Learning a Markov tree mixture).

Repeat for

i = 1, • • • , m: (a) Si = BuildMarkovTreeStructure(D) (b) θi = LearnPars(Si, D) 2. (µ) m i=1 = CompWeights((Si, θi) m i=1 , D) 3. Return µi, Si, θi m i=1

Chow-Liu algorithm

This algorithm learns a Markov tree structure maximizing the likelihood of the training set [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]). Its principle is described by Algorithm 2. It can be decomposed in two steps : Step 1. computes from the dataset the maximum likelihood estimates of the mutual informations between each pair of variables to fill an n × n symmetrical matrix (M I); Step 2. searches for a maximum weight spanning tree (MWST) in this matrix (e.g. Kruskal's algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF], used here).

Algorithm 2 (Chow-Liu algorithm).

1. M I = [0]n×n; Repeat for k = 1, • • • , n: Repeat for j = k + 1, • • • , n: i. M I[k, j] = CompM I(X k , Xj, D); ii. M I[j, k] = M I[k, j]. 2. S = CompKruskal(M I); Return S.
The first step requires O(n 2 N ) computations, while the second has a complexity of E log(E) with E the number of considered edges. In the Chow-Liu algorithm E = n(n -1)/2, so this second complexity becomes O(n 2 log(n 2 )).

Randomized edge sampling

To reduce the complexity of the Chow-Liu algorithm, we propose to apply the Perturb and Combine principle by learning each model from an incomplete matrix M I.

The random edge sampling algorithm performs this by randomly selecting a subset of a priori fixed size K of different pairs of variables according to a uniform distribution. These terms are used to partially fill the matrix M I used as input to the MWST algorithm. Algorithm 3 describes this procedure.

Algorithm 3 (randomized edge sampling).

1. M I = [0]n×n ; Repeat for k = 1, • • • , K: (a) Draw new random pair (i1, i2) ∈ {1, . . . , n} 2 ; (b) M I[i1, i2] = CompM I(Xi 1 , Xi 2 , D); (c) M I[i2, i1] = M I[i1, i2]. 2. S = CompKruskal(M I); Return S.
The complexity of Algorithm 3 is loglinear in the number K of edges drawn. Notice that the tree structure that it infers may be disconnected, and that its dependence on the dataset is increasing with the value of K. We will report in this paper simulations and results for two values of the parameter : K = n log(n) considered in [START_REF] Ammar | Subquadratic Markov tree mixture models for probability density estimation[END_REF], and K = 0.33n(n -1)/2. The first value allows a total complexity of n log(n) log(n log(n)), which is sub-quadratic and very close to the quasi-linear. The second corresponds approximately to the edges sampled by the method described in the next section.

Randomized vertex clustering

Another idea to weaken the Chow and Liu procedure was proposed by [START_REF] Schnitzler | Towards sub-quadratic learning of probability density models in the form of mixtures of trees[END_REF]. This method is a less naïve approach to sampling the matrix M I, which targets potentially interesting (i.e. of large weight) edges. Algorithm 4 details this two-step process, that first builds a local structure of the problem, and then focuses on pairs of variables located close to each other in that structure.

The first step consists in an approximate online clustering of the variables based on their mutual information, inspired by leader clustering (a cluster C p is represented by its leader L p ). As illustrated in figure 1, a sequence C of clusters is created until all variables belong to one. The construction of a cluster C p is based on two thresholds on mutual information: one cluster-threshold (M I C ) and one neighborhoodthreshold (M I N ). The cluster is built by comparing the mutual information of each remaining unclustered variable to the new leader (L p , chosen first at random, then among unclustered variables by minimizing q<p M I(L p , L q )). An unclustered variable X is identified as:

1. member of C p , if M I(X, L p ) > M I C , 2. neighbor of C p , if M I C M I(X, L p ) > M I N , 3. not related to C p , otherwise.
Setting those thresholds can be seen as excluding potentially independent variables. This exclusion rate can be controlled, since the maximum likelihood estimate of the mutual information for two independent variables asymptotically follows a χ 2 law. In this work, its percentile 0.5 (respectively 5) was used for M I C (M I N ).

Algorithm 4 (Vertex Clustering).

nbClusters = size(C);

Repeat for p = 1, • • • , nbClusters:

(a) Repeat for i1, i2 : Xi 1 , Xi 2 ∈ Cp: i. M I[i1, i2] = CompM I(Xi 1 , Xi 2 , D); ii. M I[i2, i1] = M I[i1, i2]. (b) Repeat ∀q < p : Cq ∈ N eibhors(Cp): Repeat for i1, i2 : Xi 1 ∈ Cp, Xi 2 ∈ Cq: A. M I[i1, i2] = CompM I(Xi 1 , Xi 2 , D); B. M I[i2, i1] = M I[i1, i2].
3. S = CompKruskal(M I); Return S.

In the second step of the algorithm, the mutual information of all potentially interesting pairs (X i , X j ) are computed and used as edgeweights for the MWST algorithm. Interesting pairs (a) are in the same cluster or (b) span two neighboring clusters, i.e one variable of one cluster is a neighbor of the other cluster. In addition, all edges evaluated during the clustering process are used as candidate edges.

The complexity of this algorithm is between linear and quadratic in the number of variables, depending on the numerical values of M I C and M I N and the problem structure.

Inertial search heuristic

This algorithm [START_REF] Ammar | Subquadratic Markov tree mixture models for probability density estimation[END_REF] for computing a sequence of sub-optimal MWST was designed to improve the base method from section 3.2. In this work we also apply it to the vertex clustering algorithm (section 3.3).

The inertial method takes advantage of the Markov tree structure S i-1 built in the previous iteration to partially fill the new M I i matrix. The weights (recomputed in case bootstrap copies of the dataset are used) of the edges of the Markov tree built at the previous iteration i -1 are first written in the M I i matrix of the current iteration i, and a new set of edges generated by the base method (either at random or by vertex clustering) is inserted afterwards. This is described by Algorithm 5.

The complexity of this method is similar to the base method. Algorithm 5 (Inertial research procedure).

1. M Ii = [0]n×n ; Repeat for k = 1, • • • , nbEdges(Si-1): (a) (i1, i2) = GetIndices(GetEdge(Si-1, k)); (b) M Ii[i1, i2] = CompM I(Xi 1 , Xi 2 , D); (c) M Ii[i2, i1] = M Ii[i1, i2]. 2. Repeat for k = 1, • • • , nbEdges(BaseM ethod): (a) (i1, i2) = indices of edge k from BaseMethod; (b) M Ii[i1, i2] = CompM I(Xi 1 , Xi 2 , D); (c) M Ii[i2, i1] = M Ii[i1, i2].
3. Si = CompKruskal(M Ii); Return Si.

Other variants

Two other variants (baselines) were also considered, namely random trees and bagging.

The first one draws a tree structure totally at random (i.e. independently from the dataset) through the use of Prüfer lists. Complexity is linear in n [START_REF] Ammar | High-dimensional probability density estimation with randomized ensembles of tree structured Bayesian networks[END_REF].

Bagging can also be used to increase randomization in a given method, by supplying a bootstrap replica of the original dataset to any treestructure learning algorithm. This may actually be quite productive in order to randomize tree structures (see the results below). However, as far as accuracy is concerned, it turns out to be preferable to use the full dataset for parameter estimation of each Markov tree generated by this method [START_REF] Schnitzler | Towards sub-quadratic learning of probability density models in the form of mixtures of trees[END_REF].

Empirical simulations

We apply the algorithm variants described in Section 3 to synthetic problems to assess their performance. We carried out repetitive experiments for different data-generating (or target) densities as described in Section 4.1; our results are reported in Section 4.2.

Experimental protocol

We present here results obtained on 10 different target distributions over n = 1000 binary variables, and we report them for mixtures and datasets of various sizes.

Target density generation Target densities are synthetic distribution factorizing according to a general directed acyclic graph structure. These models (structure and parameters) are generated by the algorithm described in [START_REF] Ammar | High-dimensional probability density estimation with randomized ensembles of tree structured Bayesian networks[END_REF].

Datasets We focus our analysis on rather small datasets (N = 100, 250, 1000) with respect to the number n = 1000 of variables. This replicates the usual situation in highdimensional problems, which are the motivation of this work. For each considered sample size and for each target distribution, we generate 5 different datasets.

Mixture learning For a given dataset, and for a given tree structure learning algorithm, we apply the mixture learning algorithm (Algo 1) by generating ensemble models of growing sizes (m = 1, m = 10, . . . , m = 150) in order to appraise the effect of the ensemble size on the quality of the resulting model.

In all our simulations, the parameters of the Markov tree models are learned from the dataset by maximum a posteriori estimation using uniform Dirichlet priors.

In our empirical tests we have always weighted the individual terms uniformly (i.e. µ = 1/m in Algorithm 1).

Accuracy evaluation

We assess the quality of each generated mixture by the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF], an asymmetric measure of similarity of a given distribution P T to a target distribution P G , defined by

D KL (P G || P T ) = X∈X P G (X) log 2 P G (X) P T (X)
.

Since screening all 2 1000 configurations of X is not possible, we estimate this quantity by Monte Carlo using a random sample of configurations generated according to P G :

DKL (P G || P T ) = X∼P G log 2 P G (X) P T (X)
.

In this work, we generated for each datagenerating distribution and each learning algorithm a fixed set of 50000 samples, which is then used for the Monte Carlo estimation of D KL of the models produced by the algorithm applied to the datasets issued from this distribution and with a growing number of mixture terms m.

Results and discussion

Table 1 describes the algorithm variants that we have evaluated, recalls their computational complexities, and also gives indications of their relative computing times in our implementation. The performances of these algorithms in terms of accuracy (D KL estimates) are reported in Figures 2 and3. As reference method we use the Chow and Liu single tree method (denoted by CL in the table and figures). The comparison is however not fair, because VCDU samples approximately 33% of all edges at each iteration, and ESDU only 1.4%. To provide a more accurate comparison, the latter was modified (ESDU%, in Figure 3 and Table 1) to use the same number of edges than VCDU. The results exposed in Figure 3 with this setting confirm that VCDU is also superior to ESDU%.

The inertial heuristic can be used to enhance both methods (IESDU, IVCDU) with nearly no additional complexity cost (see Table 1), leading to an increase in the improvement rate with the size of the model. Figures 2(a) and 2(b) show that IVCDU converges to a model slightly better than the CL tree (the best edges have been found, so the optimal tree is always learned), while the inertial randomized methods (IRB and IRS) with a complexity close to quasi-linear tend to approach CL when the number of mixture components grows and surpasses IVCDU which complexity is higher.

The same convergence can be observed for IESDU% (IESDU on 33% of edges instead of 1.4%) in Fig. 3, which also converges to the CL tree, albeit slower than IVCDU.

ESDU is degraded by the use of bagging (ESBU), and both methods yield performance only slightly better than random structures (MTU). We therefore conjecture that the low quality (as opposed to the low number) of edges used in ESDU introduces too much randomiza- tion to benefit from bagging. Indeed, when the quality of the edges considered increases over time (IESDU) bagging actually ameliorates the method. This hints towards that considering the use of bagging methods with CL trees or with the inertial vertex clustering could probably be improved by computing the best edges (i.e. significant on the original dataset) only once instead of repeating it for every new tree.

Experiments performed on sample sets of size 250 and 100 are reported in Fig. 2(b) and 2(c). Decreasing the number of samples from 1000 to 250 does not modify the results very much, as illustrated in Fig. 2(b). The main difference lies in the relative performance of all methods compared to CL algorithm, which seem to improve faster as m increases. This tends to indicate that randomization is increasingly more beneficial as the number of samples decreases.

The observation of the behavior of the bagging methods in Figure 2 further confirms this analysis. The gap between pairs of methods using the same algorithm on the original dataset or on bootstrap replicas (CL -CLBU, ESDU -ESBU, IESDU -IESBU) is widening when the number of samples is decreasing. Bagging of Chow-Liu trees (CLBU) is actually the best method on all dataset sizes. However, using bagging with inertial procedures and edge sam-pling yields also impressive results: the curve IESBU converges to CLBU in Fig. 2(c), while the IESBU algorithm has a much lower complexity than CLBU.

An alternative way to understand these results is in terms of over-fitting. From Fig. 2(c) we can see that IVCDU first surpasses CL after a few iterations, but the addition of subsequent terms worsens the mixture, which converges to the CL curve. Actually, at that point the model has fully learned the optimal tree, and it is repeatedly added to the model. The rise of the curve at that point signals the over-fitting. In addition, we can notice that MTU, the method using structures drawn independently from the dataset behaves better than most other methods in this context.

Observing the first tree of each model in the same figure, we can observe over-fitting again. IESBU (and ESBU) is better than IESDU (and ESDU). Likewise, the first term of CLBU is better than CL. A tree structure learned on a perturbed dataset leads to a graphical model that is more general. This behavior is still present at 250 samples, but is no longer noticeable at 1000 : IESBU and IESDU start at the same point, and the first bagged tree is worse than the optimal.

The application of BDeu weights (not reported in this paper) to the methods presented here leads to similar conclusions. But we found out that many methods actually display worse performances in terms of accuracy when combined with such Bayesian weights. This is understandable, since weighting each structure by its posterior probability makes sense for MTU only, since asymptotically only in this context the mixture will converge to a canonical Bayesian method. The methods that benefit the most from those 'Bayesian' weighting scheme are IRSBU and IRSDU. The first terms (built with few information) are gradually eliminated in favor of those identified later on (where the inertial procedure has iteratively improved the quality of the edges).

In this paper, we have compared several randomization methods aiming to approximate the Chow-Liu algorithm, with the objective of reducing its computational complexity in the context of learning mixtures of trees, and motivated by the variance reduction potential of randomization in the context of learning in highdimensional problems.

Based on our results on synthetic experiments, we claim that, in real conditions, i.e. when the number of samples is much smaller than the number of variables, randomization is interesting for probability density estimation in the form of mixtures of Markov Trees. That interest actually increases when the number of samples goes down, or when the dimensionality of the space is increasing.

In addition, we have shown that exploiting the structure of the problem by focussing on strong edges leads to methods able to compete in terms of performance with more timeconsuming procedures like bagging.

We therefore plan to keep investigating this approach. In particular, a candidate area for improvement is the transmission of knowledge between terms. Increasing the number of reused edges might speed up the convergence. Another direction of research would be the consideration of continuous variables, and the consideration of a priori known dependency/independency structures for the given problem.

  (a) First a leader (here X5) is chosen at random and compared to all 12 other variables.(b) Next, the 1st cluster is built. Here it is made of 5 members and has one neighbor.(c) The 2nd leader (X13, the farthest from X5) is compared only to 7 variables.(d) Final result, after 4 iterations. All edges considered are kept for the MWST inference.

Figure 1 :

 1 Figure 1: Illustration of the vertex clustering algorithm.

  1. V = X ; C = ∅; M I = [0]n×n; Repeat until V = ∅: (a) L = GetN ewLeader(V, C); (b) C, M I += M akeCluster(L, V, M IN , M IC ).

Figure 2 :

 2 Figure2: Average performance of the algorithms described in Table1on 5 target distributions of 1000 variables times 10 datasets, with sample sizes decreasing.

Figure 3 :

 3 Figure 3: A comparison between vertex clustering and random edge sampling methods, where the latter are modified to use the same number of edges than the first, shows that vertex clustering is still clearly superior. (1000 samples)

Table 1 :

 1 In the names of the algorithms we study, D means no alteration to the dataset and B the use of bootstrap replica. m, n, K stand for the number of terms in the mixture, of variables and of sampled pairs of variables. U emphasizes the fact that we are using uniform weights in the mixture. (CPU times are given for n = 1000 variables)

	Name	Tree generation	Dataset	Complexity	running time (one tree)
	MTU	random	D	mn	0.0017
	ESBU	rand. Edge Samp.	B	mK log(K)	0.02
	ESDU	rand. Edge Samp.	D	mK log(K)	0.02
	IESBU	Inertial Edge Samp.	B	mK log(K)	0.02
	IESDU	Inertial Edge Samp.	D	mK log(K)	0.02
	IESDU% Inertial Edge Samp.	D	(K = 165000)	0.72
	VCDU	Vert. Clust	D	up to mn 2 log(n)	0.92
	IVCDU	Inertial Vert. Clust	D	up to mn 2 log(n)	0.92
	CL	Chow-Liu	D	n 2 log(n)	1
	CLBU	Chow-Liu	B	mn 2 log(n)	1

Acknowledgments

This work was supported by FRIA/FNRS Belgium, Wallonie Bruxelles International, the French ministry of foreign and European affairs, the MESR in the framework of Hubert Curien partnerships, the BioMaGNet IUAP network of the Belgian Science Policy Office and the Pas-cal2 NOE of the EC-FP7. The scientific responsibility rests with the authors.


Sourour Ammar
email: sourour.ammar@etu.univ-nantes.fr

Philippe Leray
email: philippe.leray@univ-nantes.fr

Boris Defourny
email: boris.defourny@ulg.ac.belouiswehenkel

Density estimation with ensembles of randomized poly-trees

Motivation

Learning of Bayesian networks aims at modeling the joint density of a set of random variables from a random sample of joint observations of these variables [START_REF] Naïm | Réseaux bayésiens[END_REF]. Such a graphical model may be used for elucidating the conditional independences holding in the datagenerating distribution, for automatic reasoning under uncertainties, and for Monte-Carlo simulations. Unfortunately, currently available algorithms for Bayesian network structure learning are either restrictive in the kind of distributions they search for, or of too high computational complexity to be applicable in high dimensional spaces.

Ensembles of weakly fitted randomized models have been studied intensively and used successfully in the supervised learning literature during the last two decades. Among the advantages of these methods, let us quote the improved scalability of their learning algorithms thanks to randomization and the improved predictive accuracy the induced models thanks to their higher flexibility in terms of bias/variance trade-off. For example, ensembles of extremely randomized trees have been applied successfully in very complex highdimensional tasks, such as image and sequence classification [START_REF] Geurts | Extremely randomized trees[END_REF].

In this work we explore the Perturb and Combine idea celebrated in supervised learning in the context of probability density estimation in high-dimensional spaces. We propose a new family of unsupervised learning methods of mixtures of large ensembles of randomly generated poly-trees. The specific feature of these methods is their scalability to very large numbers of variables and training instances. We explore various variants of these methods empirically on a set of discrete test problems of growing complexity.

Methods

Poly-Tree density models

Let X = {X 1 , . . . , X n } denote a finite set of discrete random variables.

A poly-tree model P for the density over X is defined by a directed acyclic graph which skeleton is acyclic and connected, and the set of vertices of which is in bijection with X and with a set of conditional densities P P (X i |pa P (X i )), where pa P (X i ) denotes the set of variables in bijection with the parents of X i in P . It represents graphically the density factorization

P P (X 1 , . . . , X n ) = n i=1 P P (X i |pa P (X i )).
(1)

Poly-tree models can be used for probabilistic inference over P(X 1 , . . . , X n ) with a computational complexity linear in the number of variables n [START_REF] Pearl | Fusion, propagation, and structuring in belief networks[END_REF].

One can define nested subclasses of poly-tree density models by imposing constraints on the maximum number p of parents of any node. In these subclasses, not only inference but also parameter learning is of linear complexity in the number of variables. The smallest such subclass is called the tree subspace, in which nodes have exactly one parent (p = 1).

Mixture models of poly-trees

A mixture model of m poly-tree models (P 1 , . . . , P m ) is defined as a convex combination of the elementary poly-tree models, ie.

P M (X 1 , . . . , X n ) = m i=1 µ i P Pi (X 1 , . . . , X n ), (2)
where µ i ∈ [0, 1] and

m i µ i = 1.
While single poly-tree models impose restrictions on the kind of densities they can faithfully represent, mixtures of poly-trees are universal approximators.

Learning a random mixture from data

Let X be a set of discrete random variables, and

D = (x 1 , • • • , x d ) be a sample of joint observations x i = (x i 1 , • • • , x i n
) drawn from some datagenerating distribution P G (X). Let P be the space of all possible poly-tree graphical structures defined over X.

Our generic procedure for generating a random mixture of poly-tree models from D is described by Algorithm 1; it receives as inputs X, D, m, and three procedures DrawPolytree, LearnPars, ComputeWeight.

Algorithm 1 (Learning random poly-tree mixtures)

1. Repeat for i = 1, • • • , m: (a) P i = DrawPolytree(P), (b) For j = 1, • • • , n: P Pi (X j |pa Pi (X j )) = LearnPars(P i , X j , D) (c) µ i = ComputeWeight(P i , D, m) 2. Return µ i , (P Pi (X j |pa Pi (X j ))) n j=1 m i=1 .

Experiments and preliminary results

In [START_REF] Ammar | Estimation de densité par ensembles aléatoires de polyarbres[END_REF] we report some first results with the above algorithm applied to datasets of size d = 1000 generated from discrete distributions with n = 8, which could be faithfully represented by a chain, a single tree, or a single poly-tree model.

In these simulations we have considered two different instances of DrawPolytree, namely a uniform draw over the class P of all poly-trees, and a uniform draw over the subclass P 1 of trees. In order to achieve this for m ∈ {1, 2, . . . , 1000}, we have used efficient algorithms for sampling trees given in [START_REF] Quiroz | Fast random generation of binary, t-ary and other types of trees[END_REF].

For parameter learning, we used maximum a posteriori values given the dataset and structure, while assuming non-informative priors on the parameters. Concerning the µ i s, we used a uniform weighting strategy, ie.

ComputeWeight(P i , D, m) = 1/m.
Overall, these results showed that the quality of the mixture-models converges rather rapidly (ie. for m ≈ 20), and that the poly-tree mixtures were slightly superior when targeting poly-tree datagenerating distributions, while the mixtures of trees were superior in the other two cases. We also observed a slightly nonmonotonic behavior of the model quality with growing values of m, which we suspect to be related to the uniform weighting scheme.

In the immediate future, we will carry out further more systematic experiments on larger problems and spanning different versions of the algorithm.

In particular, we will consider non-uniform weighting schemes, by exploiting the score obtained for a given structure and dataset so as to downweight structures that fit less well to the datagenerating distribution. We will also consider sampling from the spaces P p of poly-trees of bounded number of parents.

Experiments will be made over a richer set of datagenerating distributions, in particular ones that can not be represented faithfully by a single poly-tree model. For instance, we will consider general directed acyclic graph models as datagenerating distributions.

We will compare our algorithm in terms of sample and computational efficiency with Bayesian network structure learning and algorithms targeting an optimal mixture of tree-models [START_REF] Meila-Predoviciu | Learning with mixtures of trees[END_REF]. Subsequently, we plan to extend our approach to handle continuous variables and incomplete datasets.

Acknowledgments

This work presents research results of the Belgian Network BIOMAGNET (Bioinformatics and Modeling: from Genomes to Networks), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office.