Annabelle Collin
email: annabelle.collin@inria.fr

Jean-Frédéric Gerbeau

Mélèze Hocini

Michel Haïssaguerre

Dominique Chapelle

J.-F Gerbeau

Surface-based electrophysiology modeling and assessment of physiological simulations in atria

come

Introduction

There is a very important medical need for modeling the electrical activity of the heart in general, and in the atria in particular, e.g. with a view to therapy planning assistance in radiofrequency ablation for patients suffering from atrial fibrillation [START_REF] Matsuo | Clinical predictors of termination and clinical outcome of catheter ablation for persistent atrial fibrillation[END_REF]. In addition to generic difficulties inherent to electrophysiology modeling, namely, modeling complexity and computational intensiveness, atria modeling features specific difficulties, in particular due to their very thin walls -mostly apparent as thick surfaces in medical imaging -which requires much refined meshes. Moreover, there is a major challenge in taking into account the anisotropy resulting from the preferred conduction direction along the muscle fibers, which is also known to vary extremely rapidly across the wall thickness [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF].

The electrical wave propagating in the cardiac tissue can be represented by a nonlinear reaction-diffusion partial differential equation (PDE), coupled with ordinary differential equations (ODEs) representing cellular activity. In this study, we consider the 2D -namely, surface-based -model proposed and mathematically substantiated in [START_REF] Chapelle | A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF], derived from the bidomain model (see for example [START_REF] Pullan | Mathematically Modeling the Electrical Activity of the Heart[END_REF][START_REF] Sachse | Computational Cardiology: Modeling of Anatomy, Electrophysiology and Mechanics[END_REF][START_REF] Sundnes | Computing the Electrical Activity in the Heart[END_REF]), and defined over the midsurface of the thin region. This surface-based model was specifically designed for thin cardiac structures -the atrial walls, in particular -and takes into account the strong anisotropy variations across the thickness. This model was already successfully numerically assessed in [START_REF] Chapelle | A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF] by comparing the resulting simulations with reference 3D simulations on thin domains of simple geometries, with dramatic benefits in computation times. Our objective here is to further assess this model with physiological simulations of the atrial electrophysiology.

Model

Atrial mesh

We produced a surface mesh representing the mid-surface of the two atria. Starting from the Zygote4 heart model -a geometric model based on actual anatomical data -we used the 3-matic5 software to obtain a computationally-correct surface mesh, and the Yams [START_REF] Frey | Yams: A fully automatic adaptive isotropic surface remeshing procedure[END_REF] meshing software to further process and refine the surface mesh. Figure 1 shows the posterior and anterior views of the mesh, which contains about 26,000 triangles and 13,500 vertices. The anatomy corresponds to ventricular end-systole, namely, when the atrial chamber has its greatest internal volume. This anatomical surface was compared with literature reports of normal human atrial dimension indicators, see for example [START_REF] Cohen | Reference values for normal adult transesophageal echocardiographic measurements[END_REF][START_REF] Jiamsripong | Three methods for evaluation of left atrial volume[END_REF], and found to be within standard ranges, see Table 1.

Crista terminalis (CT)

Bachmann's bundle (BB)

Pectinate muscles (PM)

Fossa ovalis (FO)

Right atrial isthmus

Fibers

Cardiac tissue has a fiber architecture. As the electrical conductivity is higher along than across the fiber direction, fiber orientation is very important in electric activation propagation. The specificities of the atria are that the walls are 15.0-20.0 15.9 volume (ml) 41-75 47 diameters of pulmonary veins (cm)

1.1-1.5 1.14-1.45 Table 1. Assessment of left atrial dimension indicators very thin, and that fibers orientations may vary extremely rapidly across the thickness. We use [START_REF] Ho | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF][START_REF] Krueger | Modeling atrial fiber orientation in patientspecific geometries: a semi-automatic rule-based approach[END_REF] to identify and prescribe the fibers directions at the endocardium and epicardium, see Figure 2. This figure also displays the angle θ, defined as half of the angular variation between the endocardium and the epicardium.

Surface-based bidomain model

We denote by S the midsurface of the wall. The variational bidomain surface model that we propose can be written in terms of the extracellular potential u e , the transmembrane potential V m = u i -u e , with u i the intracellular potential, as follows [START_REF] Chapelle | A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF]. Find (V m , u e ) with S u e dS = 0, such that for all t > 0,

               A m S C m ∂V m ∂t + I ion (V m ) φ dS + S σ i • ∇V m + ∇u e • ∇φ dS = A m S I app φ dS, ∀φ S (σ i + σ e ) • ∇u e • ∇ψ dS + S σ i • ∇V m • ∇ψ dS = 0, ∀ψ (1)
also with S ψ dS = 0, and where A m is a positive constant denoting the ratio of membrane area per unit volume, C m the membrane capacitance per unit surface, I (V m ) a reaction term representing the ionic current across the membrane and also depending on local ionic variables satisfying additional ODEs, and I app a given applied stimulus current.

We define the intra-and extra-cellular diffusion tensors σ i and σ e by

σ i,e = σ t i,e I + (σ l i,e -σ t i,e ) I 0 (θ)τ 0 ⊗ τ 0 + J 0 (θ)τ ⊥ 0 ⊗ τ ⊥ 0 , (2)
where I denotes the identity tensor in the tangential plane -also sometimes called the surface metric tensor -τ 0 is a unit vector parallel to the local fiber direction on the atria midsurface, and τ ⊥ 0 such that (τ 0 , τ ⊥ 0 ) gives an orthonormal basis of the tangential plane. The effect of angular variations enters the model through the coefficients I 0 (θ) = 1 2 + 1 4θ sin(2θ) and J 0 (θ) = 1 -I 0 (θ). Note that J 0 (θ) = 0 (and I 0 (θ) = 1) if and only if θ = 0 -namely, constant direction across the thickness -and then σ i,e = σ t i,e I + (σ l i,e -σ t i,e )τ 0 ⊗ τ 0 as expected for a single fiber direction. By contrast, angular variations make I 0 decrease and J 0 increase in [START_REF] Cohen | Reference values for normal adult transesophageal echocardiographic measurements[END_REF], which renders diffusion "more isotropic". This model derived from a detailed asymptotic analysis thus allows to take into account the rapid variations of the fiber direction.

The current I ion can be described by a physiological or a phenomenological model. In this study, the physiological model proposed by Courtemanche, Ramirez and Nattel in [START_REF] Courtemanche | Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach[END_REF] -most widely accepted for atria modeling -is considered.

The values of the membrane parameters are A m = 200.0 cm -1 and C m = 10 -3 mF.cm -2 . The values of the conductivity parameters vary substantially depending on the specific areas considered and are given in Table 2 (all in S.cm -1 ). The Bachmann bundle (BB), the Crista Terminalis (CT) and the pectinate muscles (PM) are regions of established fast conduction. By contrast, the Fossa Ovalis (FO) and the Isthmus of the right atrial floor (IRA) are regions of known slow conduction. All these specific regions are depicted in Figure 1.

The two atria are connected via two regions only -in the mesh connectivity, hence also in terms of electrical conduction -namely, the Bachmann bundle and the Fossa Ovalis.

Simulations results

Figure 3 displays the simulation results obtained with the above-described model.

Activation is initiated at the sinus node with a stimulus of 2 ms in duration and sufficient strength to cause the initiation of a propagating wavefront. We compare our results to several 3D modeling studies [START_REF] Krueger | Modeling atrial fiber orientation in patientspecific geometries: a semi-automatic rule-based approach[END_REF][START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF][START_REF] Deng | Simulation of biatrial conduction via different pathways during sinus rhythm with a detailed human atrial model[END_REF].

As observed in the figure, by 30 ms the wave quickly spreads along the Crista Terminalis as a consequence of the high conductivity in this part. Importantly, the depolarizing wave has now traversed the Bachmann bundle and the first activation of the left atrium has occurred at 29.6 ms. This compares very well with the findings of the modeling study [START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF] giving the left atrium activation at 29.7 ms. At t = 40 ms, we clearly see the effect of the anisotropy of the Crista Terminalis, and of the pectinate muscles already observed in [START_REF] Krueger | Modeling atrial fiber orientation in patientspecific geometries: a semi-automatic rule-based approach[END_REF], as the wavefront becomes nearly triangular -as also seen in [START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF]. The wavefront has encircled the superior Vena Cava. The Fossa Ovalis, which is the second atrial connection, becomes active (42.8 ms). By 50 ms, because of the rapid conduction in the Bachmann bundle, the wave spreads to the left atrial appendage and has activated a substantial part of the left atrial wall. The Fossa Ovalis is very active. The activation of the right atrial appendage is complete. At t = 60 ms, in the right atrium the floor is the last part that remains unexcited. The wave has encircled the mouth of the left atrial appendage at 72.8 ms. By t = 80 ms, the only part that remains unaffected by the wave in the right atrium is the isthmus, because of the reduced conductivity there. In the left atrium, at t = 90 ms we can see three separate wavefronts, also obtained in [START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF][START_REF] Deng | Simulation of biatrial conduction via different pathways during sinus rhythm with a detailed human atrial model[END_REF]. The depolarization of the left atrial appendage is complete at 95.8 ms. At 100 ms, only a small part of the left atrium in the shape of a parallelepiped is still inactive. The depolarization of the right and left atrium are complete at 101.8 ms and 109.6 ms, respectively. This is in good adequacy with the timings found in [START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF], namely, 99.3 ms and 108.2 ms, respectively, while [START_REF] Krueger | Modeling atrial fiber orientation in patientspecific geometries: a semi-automatic rule-based approach[END_REF] gives 115.0 ms for the left atrium, and [START_REF] Deng | Simulation of biatrial conduction via different pathways during sinus rhythm with a detailed human atrial model[END_REF] 119.0 ms.

For comparison purposes, we also display in Figure 4-top simulation results obtained when disregarding the fiber direction variations across the thickness, namely, when taking θ = 0. We observe major differences with the abovediscussed physiological simulation. In particular, the anisotropy effects characteristic of the Crista Terminalis are much less clearly seen, and the triangular shape is hardly observed, indeed. Overall, the activation timings are longer, in particular due to slower propagation in the Crista Terminalis, with complete depolarization occurring as late as 116.2 ms and 115.4 ms for the right and left atria, respectively. The last results displayed in Figure 4-bottom correspond to an homogeneous and isotropic tissue throughout. We increased the homogeneous conductivity with σ i = 8.0 10 -4 S.cm -1 and σ e = 3.1 10 -3 S.cm -1 to obtain roughly comparable depolarization timings. The resulting activation profile, however, is totally unrealistic -as already noted in [START_REF] Krueger | Modeling atrial fiber orientation in patientspecific geometries: a semi-automatic rule-based approach[END_REF] -see Figure 5 displaying the activation time maps for the three cases analyzed.

Concluding remarks

We have presented detailed simulation results of a surface-based electrophysiology bidomain model applied with an anatomical model of the atria. The results obtained are in excellent adequacy with previous studies, and confirm the importance of fiber-related anisotropy effects and of their strong variations across the wall thickness. Furthermore, this surface-based model provides dramatic computational benefits over 3D models with preserved accuracy [START_REF] Chapelle | A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF]. Further perspectives include applications to radio-frequency ablation planning assistance, for which computational effectiveness -less than 3 minutes for a complete simulation on a standard workstation -will be of great value.

Fig. 1 .

 1 Fig. 1. Atrial mesh posterior (left) and anterior (right) views, with main specific regions

Fig. 2 .

 2 Fig. 2. Fibers directions at the epicardium (top-left) and endocardium (top-right), angle θ (bottom-left), and angular variation visualization (bottom-right)

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Simulation of atrial depolarization

www.3dscience.com

www.materialise.com


Ahmad Ahmad Kassem

Christophe Bobineau

Christine Collet
email: christine.collet@grenoble-inp.frfuda.ma@insa-lyon.fretienne.duble

Etienne Dublé

Stéphane Grumbach
email: stephane.grumbach@inria.fr

Fuda Ma

Lourdes Martinez

Stéphane Ubéda
email: stephane.ubeda@inria.fr

Ahmad Ahmad-Kassem
email: ahmad.ahmad_kassem@inria.fr

UBIQUEST, For Rapid Prototyping of Networking Applications

Keywords: H.2 DATABASE MANAGEMENT [Languages, Systems and Software]: Query languages, Query optimisation and processing, Rule-based program execution, Distributed databases, Distributed systems, Reasoning, Information networks Algorithms, Management, Measurement, Documentation, Performance, Design, Economics, Reliability, Experimentation, Security, Human Factors, Standardization, Languages, Theory, Legal Aspects, Verification Declarative networking, programming abstraction, case-based distributed query optimization

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 trend towards ubiquitous computing is accelerated withparticularly, wireless networking -technologies interconnecting an increasing number of heterogeneous (mobile and wearable, energy constrained, personalized) devices that generate large amounts of data. These devices are autonomous, either, static or mobile and present constraints such as energy or communication capabilities. They usually take part in dedicated ad hoc networks, where applications deployment, configuration and management are tedious and require significant human involvement and expert knowledge.

In [START_REF] Ahmad-Kassem | A data-centric approach for networking applications[END_REF] we introduce our vision of a new high-level programming abstraction based on the emerging and promising declarative networking approach and declarative data manipulation expressions. Declarative networking is an emerging data-centric approach where the distributed environment is perceived as a distributed database and the applications interact through declarative queries [START_REF] Mao | On the declarativity of declarative networking[END_REF][START_REF] Loo | Declarative Routing: Extensible Routing with Declarative Queries[END_REF][START_REF] Loo | Implementing Declarative Overlays[END_REF]. This approach has been pursued at the network layer with the use of recursive query languages initially proposed to express communication network algorithms such as routing protocols [START_REF] Loo | Declarative Routing: Extensible Routing with Declarative Queries[END_REF] and declarative overlays [START_REF] Loo | Implementing Declarative Overlays[END_REF]. It has been further pursued in [START_REF] Loo | Declarative networking: language, execution and optimization[END_REF], where execution techniques for Datalog are proposed. Distributed query languages thus provide new means to express complex network problems such as node discovery [START_REF] Stillger | Leodb2's learning optimizer[END_REF], route finding, path maintenance with quality of service [START_REF] Bejerano | Algorithms for computing QoS paths with restoration[END_REF], topology discovery, including physical topology [START_REF] Bejerano | Physical topology discovery for large multi-subnet networks[END_REF], etc. The declarative networking approach is well-adapted to social systems (e.g. games, social networks, sharing), where data is pushed or pulled with incomplete knowledge in a dynamic environment.

Also declarative query languages have already been used in the context of ad-hoc networks. Several systems for sensor networks, such as TinyDB [START_REF] Madden | Tinydb: an acquisitional query processing system for sensor networks[END_REF] or Cougar [START_REF] Demers | The cougar project: a work-in-progress report[END_REF] have been proposed. They use the relational model to represent device (sensor) features and application data; they offer SQL-like languages to express data manipulation. These systems also address solutions to perform efficient data dissemination and query processing. In both cases, a distributed query execution plan is computed in a centralized manner considering the network topology and the capacities of the constrained nodes, which optimizes the placement of sub-queries in the network [START_REF] Demers | The cougar project: a work-in-progress report[END_REF][START_REF] Madden | Tinydb: an acquisitional query processing system for sensor networks[END_REF]. Declarative methods have been used also for unreliable data cleaning based on spatial and temporal characteristics of sensor data [START_REF] Jeffery | Declarative support for sensor data cleaning[END_REF].

As far as we know there is no system that integrates in a uniform way, network aspects, middleware and data management. UBIQUEST merges declarative programing languages and query languages for specifying data manipulations and distributed algorithms. Furthermore, these languages are used to intentionally express the destinations of messages, for naming and accessing data in the context of networks and dynamic environments.

The work presented in this paper describes the architecture and components of an UBIQUEST system -http://ubiquest.imag.frthat implements this approach.

It is a sort of a large distributed database system that provides a unified view of "objects" handled in networks and applications. It blurs the borders between network, operating system and middleware layers. However, from the data management point of view it should provide a means (i) to localize data (mobile applications) or define the scope of a query, (ii) to consume, filter Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. IDEAS12 2012, August 8-10, Prague [Czech Republic] Editors: Bipin C. Desai, Jaroslav Pokorny, Jorge Bernardino Copyright ©2012 ACM 978-1-4503-1234-9/12/08 $15.00. and aggregate data (continuous queries), (iii) to consider query operators that may correspond to programs, (iv) to optimize query even when no metadata or statistics are available. For that, we use Case-Based Reasoning (CBR) -learning the cost of query plans (cases) while executing them -and pseudo-random query plan generation when classical optimization techniques are inappropriate.

The paper is organized as follows. Section 2 gives an overview of our approach based on an example unifying data and networks management functions. It defines an UBIQUEST system as a set of interconnected UBIQUEST nodes. Section 3 presents the architecture of an UBIQUEST node and details its components. Section 4 focuses on the execution engines that perform program execution and global query evaluation. Section 5 presents our proof of concepts as a platform for simulating UBIQUEST systems. Finally, Section 6 concludes the paper and discusses future work.

UBIQUEST DATA-CENTRIC APPROACH

With declarative networking, the network is abstracted as a large distributed database providing unified view of "objects" handled by both networks and applications. Such a database information about the declarative programs, routers configuration, states and characteristics of the network. Rule-based programs usually correspond to network operations or protocols triggered by data updates. Rules are evaluated over local data and may communicate updates to other nodes in the network using communication primitives.

The UBIQUEST approach merges the strengths of two areas (i) databases, and (ii) declarative networking. With this approach a programmer can specify the behaviour of the system / application (the what) rather than having to describe the details of the system (the how). This allows going one step further in the overlapping approaches, for example with destinations of messages resulting of a query.

An UBIQUEST system runs on a set of computing devices interconnected through a wireless network (cf. Fig. 1). Every device embeds a virtual machine in charge of data management, processing queries (data selection and updates) and messages propagation. A message is the unit of communication among UBIQUEST nodes. It has two main parts: (i) networking information (e.g. logical destination, next hop, TTL) and (ii) a payload where the content of the message (i.e. queries or items) is embedded.

All exchanges between nodes related to communication protocols, to resource discovery or to any other applicative aspects are carried out by queries and data. This blurs the traditional distinction between communication middleware and application layers. Queries are defined using either rule-based languages (e.g. for network data query expressions or distributed algorithms) or declarative query languages (e.g. for querying application data with a global point of view). For a detailed presentation of these languages refer to [START_REF] Ahmad-Kassem | A data-centric approach for networking applications[END_REF].

Query optimization is based on CBR-based approach and pseudorandom query plan generation. This means that we learn the cost of query plans (cases) while executing them. These cases are reused for generating plans for further similar queries. If there is no convenient case, we use classical heuristics and random choice (e.g. when there is no statistics for join ordering and selection of algorithms) to generate query plans.

Figure 1. An UBIQUEST System

To illustrate our approach, let us consider an application concerning a virtual world game divided in areas and having some avatars that are located within a single area at a time (see Fig. 2). The objective of the game may be social interaction or team fighting; this does not matters for understanding the example. Every node of an UBIQUEST system has information on its own avatars and their neighbors (avatars located in the same area). Data location is thus application driven. showing that the node is the owner of the Grey avatar which is in the area 2; also the Blue avatar is in the same area but owned by node I. The Positions table at node I contains the same tuples (Fig. 2).

Such a Positions table is actually a fragment of a virtual table that maintains the information on all the avatars, giving for each of them, the area in which it is and its owner node. The virtual global Positions table for our application example is:

Considering this global view, the query to select all avatars in the virtual world and the zone where they are located is:

SELECT Area, Avatar FROM Positions;

This query can be posed at any node and the system will globally execute it. Alternatively, for restricting the scope of the query at a node level, one has to use the keyword LOCAL indicating that the query has to be evaluated over local data only.

Let us now assume that the Yellow avatar, owned by node G, is moved from area 7 (where avatar Green owned by node J is localized) to area 8 where the Red avatar (node E) is localized. The Positions table after this operation follows:

The movement is coded by several updates executed at node G (owner of Yellow) for cleaning area 7, changing the Area attributes of the avatar and finally for storing the new area exploration. The first update for cleaning the area 7 is:

Delete from Positions Where Area = (Local Select Area from Positions where Avatar = 'Yellow') and Area not in (Local Select Area from Positions where Avatar <> 'Yellow' and Owner = SELF) Stored on SELF;

The keyword LOCAL indicates that the subquery has to be evaluated by the node over local data only. The sub-queries are local and the delete operation too as it concerns only data stored on SELF. Such a query is executed at the node level and processed in a distributed way with the following principles: 1.

No centralized control. Query processing is performed in an environment that is highly dynamic, and has to adapt to and recover from the network evolution. The control needs to be fully distributed over the network. 2.

Scarce metadata. The network being highly dynamic, there is no stable knowledge on the data organization. Resource discovery is combined with networking protocols.

3.

Everything in the database. The network management is done through queries.

UBIQUEST NODE

An UBIQUEST node is a device equipped with an UBIQUEST Virtual Machine (UBIQUEST VM) complemented with a Device wrapper that allows device/VM interaction (see Fig 3). The UBIQUEST VM is composed of: (i) a Local DMS, (ii) an UBIQUEST API, and (iii) an UBIQUEST Engine comprising subengines in charge of evaluating global queries, executing rulebased programs, maintaining sensed data and the list of physical neighboring nodes.

Local DMS and UBIQUEST API

The Local DMS stores and manages data as Itemsets: application data (e.g. sensed data), network data (e.g. routing tables, neighbor table), rule-programs (e.g. distributed algorithms that can be dynamically loaded/removed to/from the system), and internal data (e.g. device specific data) used for running other UBIQUEST VM components.

The UBIQUEST API manages all interactions between the UBIQUEST Engines and the rest of the world: local applications, device sensors and other UBIQUEST VM through message exchange.

As shown in Fig. 3, the API is composed of: (I) the Application API, in charge of the interaction with applications running on the local node, (ii) the Reception and Emission modules to deal with message exchange among UBIQUEST nodes, (iii) the Sensing API that locally stores data coming from sensors embedded in the Using Payload, logical destinations and a ProgramId identifying a dissemination protocol, the Emission Module builds a new message, invokes the UNIQUEST Engine to compute the immediate physical destination(s) from the logical ones, and sends the message over the network.

The Payload Dispatcher maintains a record of the identifiers of payloads that are currently executed at the node. This allows determining if a received payload was already executed, and thus avoids loops. When it receives payloads from the Application API, it generates a new identifier for registering. When it receives payloads from the Reception Module, the Payload is forwarded to the UBIQUEST Engine for treatment.

When a payload is not in the record, the Payload Dispatcher generates a new identifier, registers the payload and transfers it to the corresponding engine. If the identifier is in the records, the payload dispatcher transfers it to the corresponding engine instance according to the query/result type and the payload identifier. If the message contains several destinations, the payload dispatcher sends it to the emission module, which constructs a message and propagates it over the network using a program (dissemination protocol) selected by the UBIQUEST Engine (Communication Module).

Sensing and Topology Engines

These two modules are autonomous and react to changes in the environment detected by the device and signaled through the Device Wrapper. The Sensing Engine gets the measures coming from physical sensors embedded in the device (e.g. temperature, location) and stores these values in corresponding itemsets. These itemsets are predefined and adopt a common structure (e.g. itemset Temperature(NodeId {key}, value)). The Topology Engine is responsible of updating the Link itemset, defined as Link(NodeId {key}, Neighbor {key}) according to physical network connections that are established or removed. The Link itemset is mandatory and is sufficient to permit communication among nodes.

Communication Module

The Communication Module has two different roles: (i) determine if the local node is part of the logical destination of incoming messages, and (ii) determine what is(are) the next hop(s) to transmit a message to a logical destination.

The logical destination of a message is either expressed extensionally using a list of node identifiers, either expressed intentionally using a query returning node identifiers, or expressed by a combination of both.

DQE Engine

The Distributed Query Engine is responsible of executing global -DLAQL queries. The DLAQL language extends the well-known SQL2 data manipulation language to conform to the data distribution policy of UBIQUEST. This means that a DLAQL expression may explicitly indicate on which UBIQUEST node data has to be stored. The Query Optimizer is based on the Case-Based Reasoning (CBR) approach as in [START_REF] Stillger | Leodb2's learning optimizer[END_REF]. It proposes to retrieve and adapt query plans using the experiences gained from the execution of past similar queries. When no knowledge is available it randomly generates query plans using classical heuristics [START_REF] Ioannidis | Query optimization[END_REF][START_REF] Subramanian | Query optimization in multidatabase systems[END_REF].

Rule Program Engines

The Rule Program Engine is in charge of executing rule-based declarative programs exploited for specifying distributed algorithms (e.g. networking protocols, sub-query execution). The engine selects which rules have to be triggered and execute them over the local data. The rule execution may involve local data storage or emission to neighborhood.

The proposed Netlog and Questlog rule-based programing languages extend Datalog with communication primitives, as well as aggregation and non-deterministic constructs which are standard in network applications. The computation of rule programs is local, and the result can be either stored locally on an UBIQUEST node on which the rules run, or sent to other nodes.

The Rule Program Engine receives payloads from the Payload Dispatcher and has to treat their Contents containing either items (new facts or query results) or predicates corresponding to queries.

If a Content contains facts, the Rule Program Engine identifies which rule-program has to be triggered by comparing new facts with predicates in the rule body (Netlog). Then, it retrieves the corresponding rules from the DMS and evaluates them in forward chaining mode till a fix point is reached.

If a Content contains a predicate corresponding to a query (Questlog), the Rule Program Engine identifies which ruleprogram has to be triggered by comparing the predicate with rule head, then it retrieves the corresponding rules from the DMS and evaluates them in backward chaining till the full query result is computed.

If a Content contains query results, these results are exploited to continue query evaluation.

Query and rule-based program execution

Distributed Query Execution

As said in section 3.4, incoming global queries are rewritten by the Scheduler as a sequence of queries that are to be evaluated in correct order to produce the final result (e.g. evaluate asynchronous subqueries before executing rewritten upper level query). These queries generally contain access to local data and to distant data, according to the horizontal fragmentation of global itemsets. An optimal query plan has to be generated for each one of these queries.

A query plan is a tree whose nodes are physical operators corresponding to data manipulation. A root node corresponds to a DLAQL command (i.e. SELECT, INSERT, DELETE, UPDATE), intermediate nodes correspond to computation operators (e.g. unions, joins, filter or aggregate), and leafs correspond to data access operators: local DMS querying, sub-query emission to all neighbors, or rule-based program invocation.

A case is generated for any (sub-)query expression Q evaluated on a node. It is composed of the expression of Q, a query plan P for Q and a set of cost parameter measures taken during the execution of P. A case is stored in the local case base of the node. The following example is a case for a query finding all avatars in area 7: {Q = "Select Avatar from Positions where Area = 7", P = Union( DMS(σ Area=7 (Positions)), SubQ(σ Area=7 (Positions))), Cost = {Energy=0.5%, Time=12ms, Memory=2%, …} } The query plan P involves computation of a subquery on local DMS (DMS operator) and emission of a global sub-query (SubQ operator). Cost parameters are normalized when possible. In our approach, the query plan generation is a recursive process.

The optimizer checks if the incoming query Q is known in the case base, i.e. if there exist similar cases corresponding to this query, select the best case/plan among them according to the optimization objectives (i.e. minimizing the given cost function). A similarity function to compare a case and a query has been defined based on classification of query expressions (i.e. based on DLAQL clauses, see [START_REF] Martinez | The QOL approach for optimizing distributed queries without complete knowledge[END_REF]). If a plan is selected then it has to be adapted to the current Q expression and to be executed, otherwise a plan is generated using a pseudo-random approach. This generation process applies classical heuristics and random choice when missing metadata (e.g. statistics) is needed. For example, if Q contains a set of join operations, it select randomly one join operation and a corresponding algorithm (join operator) and return a plan for this join operation involving two sub-queries. The same process is applied recursively for every sub-query. An in depth description of the optimization process is given in [START_REF] Martinez | The QOL approach for optimizing distributed queries without complete knowledge[END_REF].

The Execution Engine executes a query plan P using the wellknown Iterator model [START_REF] Graefe | Query evaluation techniques for large databases[END_REF] for the physical operators. It also coordinates the local and distant sub-queries and constructs a final result from sub-query results. During the execution, the cost parameters (energy, time, memory etc.) are measured and a new case is built.

Rule program execution

As we already explained the Rule Program Engine receives payloads from the Payload Dispatcher and has to treat their Contents containing either items (new facts or query results) or predicates corresponding to queries.

In addition, the Rule Program Engine propagates new items or new queries to other nodes, through the UBIQUEST API, and/or stores new items in the DMS. The Engine has some additional functions, such as timers, necessary for networking protocols, and also uses optimization techniques, such as the triggering of rules by new facts, which avoid unnecessary computations, when there are no changes in the input of rules.

Let us assume the following simple routing table maintenance protocol:

Route(SELF, dest, dest, 1) :-Link(SELF, dest). Route(SELF, dest, neigh, l2) :-Link(SELF, neigh), Route(neigh, dest, _, l1), l2 := l1 + 1. These rules will be executed if a new neighbor is discovered (i.e. new item in the Link itemset). Satisfied rules involve local storage and broadcasting of new facts (heads of the rules) to all neighbors ( symbol). This is done via the Emission Module alone as the destination is expressed extensionally (i.e. all neighbors). The following sequence diagram describes the whole process:

Combining DLAQL queries and rule programs

Going back to the evaluation of the query Q = "Select Avatar from Positions where Area = 7", as Union(DMS(σ area=7 (Positions)), SubQ(σ area=7 (Positions))); One can figure out, this leads to useless sub-query evaluation. The distant part (SubQ) is useless because the node locally stores its avatar (localized in area 7) and its neighbors (i.e. in the same area). Such knowledge can be represented as rule-based programs executing specific algorithms to solve these sub-queries. Here, this program is the local identity where no sub-query is emitted:

Positions(_, area, _) ← Positions(_, area, _).
The DQE Engine may exploit this program to solve (part of) the query Q, thanks to the correspondence table between queries and programs. The following sequence diagram shows the interactions among UBIQUEST node components for such an evaluation.

SIMULATION PLATFORM

We develop a platform (see Fig. 6) facilitating the development and monitoring of UBIQUEST applications. This platform offers tools for editing and compiling rule-based programs, a allows the simulation of network of UBIQUEST nodes.

The simulation platform has three main components:

•

The Network Editor, which allows to build and simulate a network with various UBIQUEST nodes;

•

The Network Monitor, which allows to visualize and interact with the network at run time; and •

The Node Monitor, which allows to monitor the activity of and interact with individual nodes. The Network Editor allows creating groups of nodes, displaying the status of the nodes in each group and installing rule programs on them. They can have different colors, radio range, and characteristics, such as mobile or fixed. The system creates the groups and displays the nodes on the left part of the screen. Each node is listed and for each node one can see its identifier, address, position and radio range.

The Network Monitor offers the view of the different groups of nodes, represented by different shapes and different colors, and the connections between them (if the nodes are located inside the radio range of another node). Each node has a unique identifier. The Network Monitor also allows to interact with the network, and to modify its configuration before starting or during the simulation, by moving nodes, changing their radio range, or deleting edges or nodes for instance.

The Node Monitor exhibits information about the node selected by the user, displayed on the right part of the screen. It contains six tabs: Display itemset, Programs, Messages, API, DQE, and Statistics. The Display itemset tab allows the user to choose an Itemset existing in the DMS of the node and to display the values of each items of this Itemset. It is important to notice that the content is updated on the fly. For example, you can choose to display the content of the Itemset "Route" to see all the routes contained on the selected node. The next tab simply displays which programs are installed on the node with the possibility for the user to enable or disable them on this node. We also developed a simulation and emulation environment for a detailed analysis and evaluation of queries for a large class of algorithms and protocols.

Figure 2 .

 2 Figure 2. Application scenario Such information is stored in an Itemset data structure of type: Positions(Avatar avatar{key}, Int Area, NodeID owner) For example, at node D the Position itemset has two items:

Figure 3 .

 3 Figure 3. UBIQUEST node components physical device, and (iv) the Payload Dispatcher, which manages Payload exchange among UBIQUEST VM sub-components. The Application API module validates DLAQL queries/updates submitted by applications, and translates them into an internal representation before sending them to the UBIQUEST Engine for evaluation. The Reception Module receives messages from other UBIQUEST nodes and decides if the payload of the incoming message has to be treated locally. It checks if the local node is part of the logical destination of the message. This process may involve with the UBIQUEST Engine to resolve intentional expressions of logical destinations (i.e. destination expressed using a query). Finally, the Reception Module sends the Payload of the message to the Payload Dispatcher to treat it, if the local node is one of the destinations, or forward the message to other destinations through the Emission Module, if not.

Figure 4 .

 4 Figure 4. UBIQUEST Simulation Platform displays all messages received or sent by the node. If you click on one message in particular, you can display its content. The API tab permits to modify the content of the DMS of the selected node by adding an item in one of the Itemsets of the node, or by submitting queries expressed in DLAQL, as an application would do. The DQE tab allows to monitor query execution by exploring the case base (i.e. query families, query plans and measures of computation cost), displaying the list of pending subqueries and partial results for any of the queries running on the selected node. The last tab shows some basic statistics about the node such as the number of Select queries or Update queries done in the database of the node.

  If it is expressed extensionally, determining if the local node takes part in the logical destination of a message is straightforward. In the other case, the Communication Module asks the Distributed Query Engine to solve the intentional destination (i.e. obtain extensional destinations) before deciding.

To determine the next hop(s) for propagating a message, the Communication Module selects a propagation program and invokes the Rule Program Engine to execute it. The default propagation program simply do broadcasting to all neighbors (i.e. the next hops correspond to all items of the Link itemset). Other propagation programs may be written by developers (e.g. by exploiting and maintaining a routing table) and may be automatically selected by the Communication Module.

  The Query Scheduler rewrites a global query into a set of subqueries and schedules their evaluation (e.g. a global UPDATE query is decomposed into a sequence of SELECT, DELETE and INSERT sub-queries to read the old value, delete it and insert the new value). Moreover, this module rewrites a query considering local and distant Itemset fragments generating a query (or set of queries) equivalent to the original one.

	Execution plans are composed of classical physical operators
	(implementing algebraic operations) and specific operators to
	invoke program or propagate subqueries. Efficient execution plans
	are selected using a combination of Case-Base Reasoning and
	pseudo-random query plan generation.

The role of the DQE Engine is to build and execute efficient local query execution plans according to a given cost function (expressed as a combination of real cost parameters).

The Distributed Query Engine is composed of: (i) a Query Scheduler, (ii) a Query Optimizer and (iii) an Execution Engine.

  The Messages tab

								Application API	ReceptionModule	EmissionModule	PayloadDispatcher	DQE	Rule Engine	DMS
								Query	Dispatch (Tag = DQE)	Execute Query	Execute local subquery
								Application	Results
								Emit subquery	Dispatch (Tag=Emission)
								Dispatch (Tag=DQE)	Treat results	Execute Program
								Results
								Dispatch (Tag=AppAPI)
								Results	Return results
	Application API	ReceptionModule	EmissionModule	PayloadDispatcher	DQE	Rule Engine	DMS
	Query	Dispatch (Tag = DQE)		Execute Query	Execute local subquery
	Application					Results
			Emit subquery	Dispatch (Tag=Emission)
			Dispatch (Tag=DQE)		Treat results	Execute Program
						Results
					Dispatch (Tag=AppAPI)
	Results	Return results

ACKNOWLEDGMENTS

This work has been supported by the ANR-09-BLAN-0131-01 UBIQUEST Project (http://ubiquest.imag.fr), financed by the French National Research Agency (ANR).


Pietro Gori

Olivier Colliot

Yulia Worbe

Linda Marrakchi-Kacem

Sophie Lecomte

Cyril Poupon

Andreas Hartmann

Nicholas Ayache

Stanley Durrleman

L Marrakchi

Towards joint morphometry of white matter tracts and gray matter surfaces

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

Introduction

Morphological studies of anatomical structures play an important role in neuroimaging. For example in the analysis of the circuits connecting sub-cortical areas to the cortical surface one would like to study together gray matter surfaces and white matter tracts (ber bundles). To the best of our knowledge, there is no generic method to study these objects together.

Atlas construction is based on the joint estimate of a common mean shape (template) and of its deformation towards each shape of the population. The atlas shows the invariants of the population and their variability.

We propose a new atlas construction method based on currents [START_REF] Durrleman | Statistical models on sets of curves and surfaces based on currents[END_REF] which permits to deal with curves and surfaces together. We extend to curves the methodology proposed in [START_REF] Durrleman | Topology preserving atlas construction from shape data without correspondence using sparse parameters[END_REF] for surfaces. This allows to x a topologically correct representation of the bundle templates. We are therefore able to compare the templates directly with the shapes of the population and to study the relative positions of dierent objects.

Moreover every deformation is based on one single dieomorphism of the whole 3D space which preserves the spatial organization of the objects making possibile a joint analysis of multiple objects.

Methods

Given a population of N subjects, the aim of this method is to estimate simultaneously the template complex and its deformations towards the N shape complexes. These complexes gather together all the objects under examination.

The number of deformations is equal to the number of subjects. Every deformation is based on one single dieomorphism which is built using a nite set of control points shared among the population. The deformation parameters are given by initial momenta linked to the control points. Thus, the number of control points denes the dimension of the parameterization of the deformations.

The method proposed here is based on the minimization of a cost function achieved with a gradient descent method. We need therefore an initial template complex. We choose an ellipsoid for the surfaces and we propose a new generic method for the bundle templates by selecting the most representative bers among all the subjects present in the population with a greedy approximation method based on the framework of currents (see Fig. 1).

This diers from a previous atlas method for curves based on currents [START_REF] Durrleman | Registration, atlas estimation and variability analysis of white matter ber bundles modeled as currents[END_REF] where the templates were not given as ber bundles, thus making dicult the analysis of the relative position between dierent objects.

We apply our new method on two populations, one of 5 controls and one of 5 patients with Gilles de la Tourette syndrome [START_REF] Worbe | Distinct structural changes underpin clinical phenotypes in patients with gilles de la tourette syndrome[END_REF]. For each population the data set consists of the left globus pallidus (lGP) and the bundle of bers connecting the globus pallidus to the cortical surface (BGP). Their segmentation was performed from T1-weighted MRI and DWI respectively [START_REF] Worbe | Altered structure of cortico-striato-pallido-thalamic circuits in gilles de la tourette syndrome[END_REF].

Results

The nal estimated template complex for both populations and the initial template complex are shown in Fig. 1. The topology of the initial template complex has been preserved throughout the atlas construction. The comparison between the two updated template complexes reveals the dierences between the common features of the two populations.

In Fig. 2 we show the updated template complex of the population of controls with the initial momenta of each subject.

We perform also two PCA on the deformation parameters of each population. Fig. 3 shows the rst mode for both populations at ±σ. The main variability in the controls is represented by a torque of BGP and lGP towards the posterior part of the brain while in the patients there is mainly a shrinkage/elongation of both BGP and lGP in the upper part of the lGP.

Conclusions

We proposed a new atlas construction method which permits the joint analysis of white matter tracts and gray matter surfaces in a way which preserves their topology in the templates. This method allows the study of the relative positions of dierent objects while respecting their underlying anatomical organization.

Figure 1

 1 Figure 1

Figure

  Figure 2

2

 2 Figure 2

3

 3

Figure 3

 3 Figure 3